APPLICATION FOR UNITED STATES LETTERS PATENT 

f 



INVENTORS: Weidong ZHU, Guangyao XU, Chun-Nam WONG, and Nengan 
ZHENG 



TITLE: SYSTEM AND METHOD FOR DETECTING STRUCTURAL 

DAMAGE 



ATTORNEYS: 
& 

ADDRESS: 



FLESHNER & KIM, LLP 
P. O. Box 221200 
Chantilly,VA 20153-1200 



DOCKET NO.: UMBC-0015 




SYSTEM AND METHOD FOR DETECTING STRUCTURAL DAMAGE 

BACKGROUND OF THE INVENTION 

This application claims priority to U.S. Provisional Patent Application Serial No. 
60/471,873 filed May 20, 2003 and U.S. Provisional Patent Application Serial No. 
60/512,656 filed October 10, 2003, which are both incorporated herein by reference. 

1. Field of the Invention 

[1] The invention relates to a method and apparatus for detecting structural 
damage, and, more specifically, to a method and apparatus for detecting structural damage 
using changes in natural frequencies and/ or mode shapes. 

2. Background of the Related Art 

[2] Damage in a structure can be defined as a reduction in the structure's load 
bearing capability, which may result from a deterioration of the structure's components and 
connections. All load bearing structures continuously accumulate structural damage, and 
early detection, assessment and monitoring of this structural damage and appropriate 
removal from service is the key to avoiding catastrophic failures, which may otherwise result 
in extensive property damage and cost. 

[3] A number of conventional non-destructive test (NDT) methods are used to 
inspect load bearing structures. Visual inspection of structural members is often 
unquantifiable and unreliable, especially in instances where access to damaged areas may be 
impeded or damage may be concealed by paint, rust, or other coverings. Penetrant testing 
(PT) requires that an entire surface of the structure be covered with a dye solution, and then 
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inspected. PT reveals only surface cracks and imperfections, and can require a large amount 
of potentially hazardous dye be applied and disposed of. Similarly, magnetic particle testing 
(MT) requires that an entire surface of the structure be treated, can be applied only to 
ferrous materials, and detects only relatively shallow cracks. Further, due to the current 
required to generate a strong enough magnetic field to detect cracks, MT is not practically 
applied to large structures. Likewise, eddy current testing (ET) uses changes in the flow of 
eddy currents to detect flaws, and only works on materials that are electrically conductive. 
Ultrasonic testing (UT) uses transmission of high frequency sound waves into a material to 
detect imperfections. Results generated by all of these methods can be skewed due to 
surface conditions, and cannot easily isolate damage at joints and boundaries of the 
structure. Unless a general vicinity of a damage location is known prior to inspection, none 
of these methods are easily or practically applied to large structures which are already in 
place and operating. On the other hand, resonant inspection methods are not capable of 
determining the extent or location of damage, and is used only on a component rather than a 
assembled structure. None of the above NDT methods are easily or practically applied to 
large structures requiring a high degree of structural integrity. 

[4] Because of these shortfalls in existing NDT methods when inspecting 
relatively large structures, structural damage detection using changes in vibration 
characteristics has received much attention in recent years. Vibration based health 
monitoring for rotating machinery is a relatively mature technology, using a non-model 
based approach to provide a qualitative comparison of current data to historical data. 
However, this type of vibration based damage testing does not work for most structures. 
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Rather, vibration based damage detection for structures is model based, comparing test data 
to analytical data from finite element models to detect the location(s) and extent of damage. 
Vibration based damage detection methods fall into three basic categories. The first of these 
is direct methods such as optimal matrix updating algorithms, which identify damage 
location and extent in a single iteration. Because of the single iteration, these methods are 
not accurate in detecting a large level of damage. The second category is iterative methods. 
The methodology has only been for updating modeling, which determines modified 
structural parameters iteratively by minimizing differences between model and test data. The 
third category includes control-based eigenstructure assignment methods, which have the 
similar limitation to that of the direct methods indicated above and are not accurate in 
detecting a large level of damage. None of these current vibration based methods have been 
incorporated into an iterative algorithm that can detect small to large levels of damage, and 
the vibration based approach for structures remains an immature technology area which is 
not readily available on a commercial basis. 

SUMMARY OF THE INVENTION 

[5] An object of the invention is to solve at least the above problems and/or 
disadvantages and to provide at least the advantages described hereinafter. 

[6] Another object of the invention is to provide a system and method for 
detecting structural damage based on changes in natural frequencies and/or mode shapes. 
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[7] An advantage of the system and method as embodied and broadly described 
herein is that it can be applied to a large operating structure with a large number (thousands 
or more) of degrees of freedom. 

[8] Another advantage of the system and method as embodied and broadly 
described herein is that it can accurately detect the location(s) and extent of small to large 
levels of damage and is especially useful for detecting a large level of damage with severe 
mismatch between the eigenparameters of the damaged and undamaged structures. 

[9] Another advantage of the system and method as embodied and broadly 
described herein is that it can work with a limited number of measured vibration modes. 

[10] Another advantage of the system and method as embodied and broadly 
described herein is that it can use measurement at only a small number of locations 
compared to the degrees of freedom of the system. A modified eigenvector expansion 
method is used to deal with the incomplete eigenvector measurement problem arising from 
experimental measurement of a lesser number of degrees of freedom than that of the 
appropriate analytical model. 

[11] Another advantage of the system and method as embodied and broadly 
described herein is that it can be applied to structures with slight nonlinearities such as 
opening and closing cracks. The random shaker test or the random impact series method 
can be used to average out slight nonlinearities and extract linearized natural frequencies 
and/ or mode shapes of a structure. 
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[12] Another advantage of the system and method as embodied and broadly 
described herein is that it can handle structures with closely spaced vibration modes, where 
mode switching can occur in the damage detection process. 

[13] Another advantage of the system and method as embodied and broadly 
described herein is that it can handle different levels of measurement noise with estimation 
errors within the noise levels. 

[14] Another advantage of the system and method as embodied and broadly 
described herein is that the damage detection method and the vibration testing methods 
such as the random impact series method enables damage detection and assessment to be 
automated, thus improving the reliability/integrity of results. 

[15] Another advantage of the system and method as embodied and broadly 
described herein is that damage detection and assessment may be automated in the field so 
that structural health can be monitored at central location and useful service life may be 
optimized. 

[16] Another advantage of the system and method as embodied and broadly 
described herein is that the random impact series method enables the modal parameters such 
as natural frequencies and/or mode shapes to be measured for a large structure or a 
structure in the field when there are noise effects such as those arising from the wind or 
other ambient excitation. 

[17] Additional advantages, objects, and features of the invention will be set forth 
in part in the description which follows and in part will become apparent to those having 
ordinary skill in the art upon examination of the following or may be learned from practice 
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of the invention. The objects and advantages of the invention may be realized and attained 
as particularly pointed out in the appended claims. 

BRIEF DESCRIPTION OF THE DRAWINGS 

[18] The invention will be described in detail with reference to the following 
drawings in which like reference numerals refer to like elements wherein: 

[19] Figure 1A shows a system 100 for detection of structural damage according to 
one embodiment of the invention; 

[20] Figure IB is a flowchart of an inverse algorithm for identifying stiffness 
parameters of a damaged structure from a select set of measured eigenparameters, in 
accordance with an embodiment of the invention; 

[21] Figure 2 is a flowchart of a quasi-Newton method for finding an optimal 
solution to the system equations shown in the flowchart of Figure IB, in accordance with an 
embodiment of the invention; 

[22] Figure 3 is a schematic view of a serial mass-spring system; 

[23] Figures 4A-4C illustrate estimation errors in a series of iterations for a low 
order system with a small level of damage; 

[24] Figures 5A-5C illustrate estimation errors in a series of iterations for a low 
order system with a large level of damage; 

[25] Figures 6A-6C illustrate estimation errors in a series of iterations for a large 
order system with a large level of damage; 

[26] Figure 7 illustrates a finite element model of a fixed-fixed beam; 
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[27] Figures 8A-8C illustrate the estimated stiffness parameters with complete 
eignenvector measurements and different noise levels for a ten element beam with a large 
level of damage; 

[28] Figures 9A-9B illustrate stiffness parameters with reduced eigenvector 
measurements and different noise levels for a ten element beam with a large level of damage; 
[29] Figure 10 illustrates a modular, four bay space frame; 

[30] Figure 11 illustrates the estimated dimensionless stiffnesses for the damaged 
frame as shown in Figure 11; 

[31] Figure 12 illustrates a cantilever aluminum beam test specimen with uniform 
damage in approximately 5 elements; 

[32] Figure 13 illustrates the experimental results for the estimated bending 
stiffnesses of all the elements of a cantilever aluminum beam test specimen with the actual 
damage between 10 cm and 15 cm from the cantilevered end using a 40-element finite 
element model; 

[33] Figure 14 illustrates the experimental results for the estimated bending 
stiffnesses of all the elements of a cantilevered aluminum beam test specimen with the actual 
damage between 25 cm and 30 cm from the cantilevered end using a 40-element finite 
element model; 

[34] Figure 15 illustrates the experimental results for the estimated bending 
stiffnesses of all the elements of an undamaged cantilever aluminum beam test specimen 
using a 40-element finite element model; 
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[35] Figure 16 illustrates a cantilever aluminum beam test specimen with a narrow 

cut; 

[36] Figure 17 illustrates the experimental results for the estimated bending 
stiffnesses of all the elements of the test specimen as shown in Figure 16 using a 45-element 
finite element model; 

[37] Figure 18 illustrates the estimated bending stiffnesses of all the elements of the 
cantilever aluminum with simulated damage between 9 cm and 15.75 cm from the 
cantilevered end using a 40-element finite element model; 

[38] Figure 19 illustrates the estimated bending stiffnesses of all the elements of the 
cantilever aluminum beam with simulated damage between 29.25 cm and 36 cm from the 
cantilevered end using a 40-element finite element model; 

[39] Figure 20 illustrates the estimated bending stiffnesses of all the elements of the 
cantilever aluminum beam with multiple simulated damage — one at the 3 rd element and the 
other at the 20 th element; 

[40] Figure 21 illustrates the estimated bending stiffnesses of all the elements of the 
cantilever aluminum beam with simulated uniform damage from the 16 th to the 18 th element; 

[41] Figure 22 illustrates a 50-foot lightning mast test specimen with an eccentric 
spike at the top; 

[42] Figure 23 illustrates the experimental results for the estimated bending 
stiffnesses of all the elements of the lightning mast as shown in Figure 22; 
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[43] Figure 24 illustrates the estimated bending stiffnesses of all the elements of the 
lightning mast as shown in Figure 22 with simulated damage between 0.76 m and 1.125 m 
from the ground using a 40-element finite element model; 

[44] Figure 25 illustrates the estimated bending stiffnesses of all the elements of the 
lightning mast as shown in Figure 42 with simulated damage between 6.89 m and 7.35 m 
from the ground using a 40-element finite element model; 

[45] Figures 26A-26C illustrate the estimated bending stiffnesses of all the elements 
of a cantilever aluminum beam using the first regularization method; 

[46] Figures 27A-27C illustrate the estimated bending stiffnesses of all the elements 
of a cantilever aluminum beam using the second regularization method; 

[47] Figure 28 illustrates an example of a practical application of the damage 
detection method of Figure 1A, in accordance with an embodiment of the invention; 

[48] Figure 29 illustrates the application of a series of random impacts to the 
practical application of the damage detection method of Figure 1A, in accordance with an 
embodiment of the invention; 

[49] Figure 30 illustrates a system for applying the series of random impacts shown 
in Figure 29 to the practical application of the damage detection method of Figure 1A, in 
accordance with an embodiment of the invention; 

[50] Figure 31 is a block diagram of one preferred embodiment of the random 
impact device of Figure 30; 

[51] Figure 32 illustrates a lightning mast specimen; 
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[52] Figures 33A-33B are graphical representations of coherence and FRF from 
single and multiple impact tests of the mast shown in Figure 32; 

[53] Figure 34 illustrates a random series of pulses with the same deterministic 
shape and random amplitudes and arrival times; 

[54] Figure 35 illustrates an average normalized shape function of force pulses; 

[55] Figures 36-38 are graphical representations of analytical and numerical 
solutions; 

[56] Figure 39 shows the comparison of the probability density function (PDF) of 
the experimental arrival time with that from uniform distribution; 

[57] Figure 40 shows the comparison of the experimental and analytical PDFs for 
the number of arrived pulses; and 

[58] Figure 41 shows the comparison of the experimental and analytical PDFs for 
the pulse amplitudes. 

DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS 

[59] Commonly measured modal parameters, such as natural frequencies and mode 
shapes, are functions of physical properties of a particular structure. Therefore, changes in 
these physical properties, such as reductions in stiffness resulting from the onset of cracks or 
a loosening of a connection, will cause detectable changes in these modal parameters. Thus, 
if the changes in these parameters are indicators of damage, vibration based damage 
detection may be, simplistically, reduced to a system identification problem. However, a 
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number of factors have made vibration based damage detection difficult to implement in 
practice in the past. 

[60] The system and method for detecting structural damage as embodied and 
broadly described herein is motivated by the observed advantages of vibration based damage 
detection over currently available technologies. It is well understood that this system and 
method may be effectively applied to damage detection and assessment for substantially all 
types and configurations of structures, including, but not limited to, simple beams, hollow 
tubes, trusses, frames, and the like. However, simply for ease of discussion, the system and 
method will first be discussed with respect to three examples - a mass-spring model, a beam, 
and a space frame - for conceptuali2ation purposes. The system and method will later be 
applied lightning masts in electric substations. 

[61] Figure 1A shows a system 100 for detecting structural damage according to 
one embodiment of the invention. System 100 includes a stiffness parameter unit 103 for 
detecting stiffness of a structure in question. Stiffness parameter unit 103 may include a 
spectrum analyzer or any other device capable of receiving vibration information and 
providing natural frequency and/or mode shape data. One example might be a spectrum 
analyzer capable of performing a fast fourier transform (FFT) and with modal analysis 
software for deriving mode shape data from vibration information received from sensor 110 
if mode shape information is needed. A roving sensor technique or multiple sensors may 
need to be used if mode shape information needs to be used. If mode shape information is 
desired, a sensor should be provided at the tip of the hammer or device used to excite the 
structure. This sensor should be configured to send data to the spectrum analyzer in order to 
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obtain the frequency response function, which the spectrum analy2er (modal analysis 
software) uses to obtain mode shape information. Stiffness parameter unit 103 is configured 
to receive vibration information output from the sensor 110 via sensor coupler 113 and 
input 114. Sensor coupler 113 could be a wire, optical fiber or even a wireless connection 
between sensor 110 and stiffness parameter unit 103. 

[62] Stiffness parameter unit 103 may include an iterative processing unit 115 
capable of determining stiffness parameters using a first order perturbation approach and the 
generalized inverse method. The first order perturbation approach can include using natural 
frequencies and/or mode shapes of the structure according to a preferred embodiment of 
the invention. Iterative processing unit 115 may be capable of converging to a correct result 
for the output of stiffness parameters even when the system is underdetermined. Stiffness 
parameter unit 103 may further include an outer iterative processing unit 117 and an inner 
(nested) iterative processing unit 119 which may operate using a first or higher order 
perturbation approach and a gradient or quasi-Newton method to be discussed below. 
Stiffness parameter unit 103 outputs stiffness parameters to damage information processor 
123. Damage information processor includes a damage extent processor 125 and a damage 
location processor 127. Damage location processor 127 outputs the location(s) of damage 
and damage extent processor 125 outputs the extent of damage. 

[63] System 100 using the iterative processing units 115 or 117 and 119 is capable 
of determining the stiffness parameters and ultimately the location(s) and extent of damage 
for structures with a largest dimension less than 50 meters, 25 meters, 10 meters, 2.5 meters 
and even 1.5 meters. Examples of this are provided herein. 
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[64] The system and method may include a multiple-parameter, general order 
perturbation method, in which the changes in the stiffness parameters are used as the 
perturbation parameters. By equating the coefficients of like-order terms involving the same 
perturbation parameters in the normalization relations of eigenvectors and the eigenvalue 
problem, the perturbation problem solutions of all orders may be derived, and the 
sensitivities of all eigenparameters may be obtained. The perturbation method may be used 
in an iterative manner with an optimization method to identify the stiffness parameters of 
structures. 

Methodology 

[65] This method presented below can simultaneously identify all the unknown 
stiffness parameters and is formulated as a damage detection problem. Since the effects of 
the changes in the inertial properties of a damaged structure are usually relatively small, only 
the changes in the stiffness properties due to structural damage are considered. 

[66] Consider a A^-degree-of-freedom, linear, time-invariant, self-adjoint system 
with distinct eigenvalues. The stiffness parameters of the undamaged structure are denoted 

by G w (/ = 1,2 m), where m is the number of the stiffness parameters. Structural damage is 

characterized by reductions in the stiffness parameters. The estimated stiffness parameters 
of the damaged structure before each iteration are denoted by G, (i=l,2,...,m), and its 
stiffness matrix, which depends linearly on G,., is denoted by K = K(G) , where 
G=[G l ,G 2 ,--,G m ] T . Here the superscript T denotes matrix transpose. The eigenvalue 
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problem of the structure with stiffness parameters G i is 

K# k =A, k M0 k (1) 
where M is the constant mass matrix, and X k =A*(G) and $ k =<p k (G) (k = 1,2,...,N) are the k- 

th eigenvalue and mass-normalized eigenvector respectively. It is noted that A k -co\^ where 
co k is the £-th natural frequency of the structure. The normalized eigenvectors of (1) satisfy 
the orthonormality relations 

tffMf^S^ (t k ) T Kf=A k S ku (2) 

where \ <u<N and is the Kronecker delta. Before the first iteration, the initial stiffness 

parameters of the damaged structure are assumed to be G- 0) =<x i G jW (i = l,2,...,m), where 
0 <(j. <1, and the eigenvalue problem in (1) corresponds to that of the structure with 
stiffness parameters G\ 0) . If there is no prior knowledge of the integrity of the structure, 
one can start the iteration from the stiffness parameters of the undamaged structure and set 
a i =1. Let G di (/ = l,2,...,m) denote the stiffness parameters of the damaged structure. The 
eigenvalue problem of the damaged structure is 

K/=« (3) 
where K f/ =K(G,,) is the stiffness matrix with = [G d{ ,G d2 ,~-,G dm f y and A d =A k (G d ) and 
$ d = <j> k {G d ) are the k-th eigenvalue and mass-normalized eigenvector respectively. The 
stiffness matrix K rf is related to K through the Taylor expansion 

K d =K(G,) = K + £— SG, (4) 

i=l OUi 
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where SG i =G (li ~G i (i = 1,2,..., aw) are the changes in the stiffness parameters, and the higher- 
order derivatives of K with respect to G i vanish because K is assumed to be a linear 
function of G r Based on the finite element model, the global stiffness matrix of a 
distributed structure satisfies (4) as its higher-order derivatives with respect to each element 
stiffness parameter vanish. 

[67] Let the k-xh eigenvalue and mass-normalized eigenvector of the damaged 
structure be related to A k and <f> k through 

m mm m m m 

* + Z V G < +11^*^ + ^.JG.SGy-SG, + el (5) 

''=1 '=1 7=1 /=1 j=\ /=1 

v v < 

p summations 

m m m m m m 

^=^+I<^+Zi:<,^G 7+ ... + XZ---S *l )l ,..,SG i SG,..SG l+ e> (6) . 

p summations 

where ^* y , ^* 2) .., and Af p) .,^ are the coefficients of the first, second, and p-th order 
perturbations for the eigenvalue, z* )M z k (2)ij , and z*^.,., are the coefficient vectors of the 
first, second, and p-th order perturbations for the eigenvector, and e\ and e* are the 
residuals of order p+\. Note that the numbers in the parentheses in the subscripts of the 
coefficients and coefficient vectors indicate the orders of the terms. By the Taylor 
expansion, one has for any p > 1 , 
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By (7), and z ( * p)ff . w are symmetric in the indices, i 9 j\ . .., and /. The right-hand sides 

of (7) are the ^>-th order sensitivities of the eigenvalues and eigenvectors with respect to the 
stiffness parameters. 

[68] Using the normalization relations of the eigenvectors, </> k and and 
symmetry of the coefficient vectors in (6), as indicated earlier, one obtains 



i=i i=i j=i i=i j=\ r=i 



/> summations 
mm m 



1=1 >=1 /=1 7=1 f=l 



p summations 



m m m 1 

= i + 2wvmz; 1);+ (z; 1v )W]^^ 

wi m m i 

+ 2!(z; 2) , ; ) r M^}^,^ +III^{3WYmz; 3) , ( + 2![(z; 1)j ) r Mz; 2)y , ♦(^/mz^, 

'=1 7=/ /-/ K ijl 

+ (^) r Mzl )u ) + 2\[(^ 

_1_ 

.i 

■>•■ 



+ - + ZS-S 7r{p!^) r %>, +(/'-0![(z; i) ,) r Mz; /) . 1)y/ ...„ + (z} )> ) r Mzf p . 1)tf ... J , 

'"=1 J=/ A /> -5/ 



p summations 

+ • • • + (z*,, / Mz^,.,, ] + 2 !(/> - 2) ! [(zf 2) .. f M Z ; H) ,, fl + (zf 2)j/ ) T Mz k {p _ 2) ... si 



+ • • • + (zf 2) Mz^. . J + • • • + ip - 2) ! 2 ! [(zf^,..,, ) T Mz\^H^ p . 2) ,.J^U, 



+ - + (<,-,),.., ) r Mz;, )( ] + p !(z; p) ,.. w / MjJ* }<?G,*G, • -60,60, + • ■ • (8) 
where the superscript T denotes transpose. For any />-th ( /> > 1 ) order term in the last 
expression of (8), the coefficient R*..., is defined by R*..., = A r ,!A r 2 !---A' a !, where A',, X 2 , ...... , 

and Z 0 (l<o<p) are the numbers of the first, second, and last distinct indices within 
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indices, ij, and /, with X, + X 2 + - + X a = p . For instance, fl 1234 = 1!1!1!1! = 1 with a = 4, 
and ^i 122 23 = 2 !3!1! = 2!3! with a = 3. Some explanations of the general p-th order term in the 
last expansion in (8) are in order. It includes p+\ types of terms ordered from the beginning 
to the end of the expression within the braces, with each type of terms except the first and 
last ones enclosed in square brackets. The c-xh {\<c< p + \) type of terms is obtained by 
multiplying a (M)-th order term in the expansion of ($) T by a (£-rH)-th order term in the 
expansion of M$ . For the *-th type of term, where 2 < c < p , the p indices, /,yV", and t 
are distributed in the subscripts of the two vectors pre- and post-multiplying M, whose 
numbers of indices in the subscripts are cA and p-c+1 respectively. For the first and last 
((/>+l)-th) types of terms, all the p indices lie in the subscripts of the vectors post- and pre- 
multiplying M , respectively. The number of terms within each set of square brackets equals 
the number of different combinations of indices in the subscripts of the vectors pre- and 
post-multiplying M. When all the p indices, i 9 j 9 — 9 and t y have distinct values, due to 
symmetry of these vectors in their indices, terms of the <:-th (l<c$/? + l) type, involving 
different permutations of the same indices in the subscripts of the vectors, are equal and 
combined, resulting in the scaling factor (c-l)!(/?-c + l)! in front of the square brackets. 
Consequently, the indices in the second through ^>-th summations range from the values of 
their preceding summation indices to m. When any of the p indices, i 9 j,-~, and t, have 
equal values, the corresponding terms in each type are given by those in the previous case 
divided by 7^...,. For instance, the 4 th order term of the form SG^G^ in the expansion of 
{(/> d k ) T M(f> d k can be obtained from the expression for the p-th order term in (8): 
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y^WYMz^ + 3![(z; 0l ) 7 Mz; 3)222 + (z5 )2 ) r Mzf 3)122 + (zj )2 )^; 3)122 + (zJ )2 ) r Mzf 3)122 ] 

+ 2!2![(z; 2)12 ) r M Z ; 2)22 +(z; 2)12 ) 7 Mzf 2)22 +(z; 2)12 fMz; 2)22 +(z; 2)22 ) r Mz; 2)12 +(z; 2)22 ) r Mz; 2)12 

^-(z^fMz^J + SIKz^fMzJj, + (z; 3)122 ) r Mzf 1)2 + (zf 3)(22 ) r Mz} )2 + (zf 3)122 ) r Mz} )2 ] (9) 
+ 4!(zf 4)1222 ) 7 -M^^ 

+ 4(zf 2)12 ) 7 'Mzf 2)22 + [(z; 3)222 ) r Mz; 0l + 3(zf 3)122 ) 7 'Mz; i)2 ] + 4(z; 4)1222 ) 7 M^ i }^G,^ 
where p = 4 and the four indices involved, i 9 j 9 /,and o , are / = 1 and 7 = / = 0 = 2 . 

[69] Equating the coefficients of the SG ( (1 = 1,2,..., in) terms in (8) and using 
symmetry of the mass matrix yields 

(#*) r Mkj y =0 (10) 
Equating the coefficients of the SGfiGj terms and using symmetry of M and z* 2) .. yields 

(^Mz^-icz'/Mzf,), (11) 
for all ij = l,2,...,w . Following a similar procedure, one obtains 

u k ) T Mz k {3)iJI =-|[ (z;, ) ,.) r Mz; 2) . /+ (z* ) .) 7 -Mz; 2); ,. + (z; 1)/ ) r Mz[ 2) ,] ( i 2 ) 

for i,j,l = l,2,...,m. Equating the coefficients of the SGfiGj-SG^G, terms ofp-th order in 
(8) yield 

(#*rM^--^{(p-i)i^ 

+ 2! ( ^-2)![(z; 2) ..) 7 mz; p . 2) ,.., ( +(^ r Mzf,_ 2) ,...^ (13) 

+ (/> - 2)!2! [(z^,.,, ) 7 ' Mi* „ + (z^,...,, ) T Mzf 2)j/ + - + (z^.,.., ) 7 ' Mz\ 2>l ] 

for i,j,---,t = l,2,...,m. The — 1 types of terms, enclosed in the p-1 sets of square brackets 
on the right-hand side of (13), are ordered from the beginning to the end of the expression 
within the braces, and their structures are readily observed. When p is odd, by symmetry of 
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Mand z* p) . w , thej-th (1 < y < ^— ) type of terms is identical to the (p-y)-th type, and the 

two types of terms can be combined. Similarly, when p is even, the j/-th (1 < y < y - 1 ) type 
of terms equals and can be combined with the (p-y)-th type of terms. In this case, 
however, there is a separate type, the (^) -th type, of terms in the middle of the expression. 
[70] Substituting (4)-(6) into (3) yields 



m m m 



m pjVf m m m 

* + Z^sg, ? + z«J/c, + ZZ<** W + Z Z Z z o)iji SG i 5G j SG i +• 



- * +Z^ +YLKvj sg < sg j +ZZZlV G W G ' +••• M < 14 ) 

[ /=! /=1 7=1 /=1 y=l /=! J 

f"i mm m m m ] 

*" +Z<^, +ZS"(VW +ZZZ«o^W +■■■ 
/=l 1=1 y=l /=! y=l M J 

Equating the coefficients of the <SG,. (/ = 1,2, ...,/«) terms in (14) yields 

Kl o> + ^* =A * Mz o> + 4 M ^ ( 15 ) 

Expanding z*^. using normalized eigenvectors of (1) as basis vectors, one has 

4> = ZC*" (16) 

u=l 

where P* )iu is the coefficient of the w-th term and the number in the parentheses in its 
subscript corresponds to that of zf 1)r Pre-multiplying (15) by (# k ) T and using (1), (2), and 
(16) yields 

4 -<rr±f (it) 

Substituting (16) into (10) and using (2) yields 
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/>',,=() (18) 



Pre-multiplying (15) by (f) T , where l<v<W andv* * , and using (1), (2), and (16) yields 

p* = 1 (f) T —f (19) 

By (16), (18) and (19) we have determined z{ 1)r 

[71] Equating the coefficients of the 8G i 8G } terms in (14) yields 

2!Kz; 2) , + Hz;,, + f^, =2U<Mzf 2 , + «, + + 2!^M#* (20) 

Expanding z* 2) .. using normalized eigenvectors of (1) as basis vectors, one has 

where /f 2)ffa is the coefficient of the w-th term and the number in the parentheses in its 
subscript corresponds to that of zj 2) ... Pre-multiplying (20) by (<j) k ) T and using (1), (2), (10), 
and (21) yields 

*.->«f to (22) 

Substituting (21) into (11) and using (2) yields 

$*-far"4>, (23) 

Pre-multiplying (20) by (f) T , where 1 < v < N and v * k , and using (1), (2), and (21) yields 
By (21), (23) and (24) we have determined z[ 2) ... 
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[72] Equating the coefficients of the SG^GjSG, terms in (14) yields 

VKJ i UKJj ukj { ^25) 

Expanding z* 3)iy/ using normali2ed eigenvectors of (1) as basis vectors, one has 

n 

4)i/' = Z^W" (26) 

where /^ } . 7w is the coefficient of the w-th term and the number in the parentheses in its 

subscript corresponds to that of z[ 3) . 7 . Pre-multiplying (25) by (# k ) T and using (1), (2), (10), 

and (26) yields 

2! 3K 5K 9K 

*v* *Ti i ' kn eG l t ™ +~dG*w + dGi z ™ "^(V -^(V -4 Mz (%] ( 27 ) 

Substituting (26) into (12) and using (2) yields 

J&» =-|}[ ( z w) r ^(V + ^)y) r ^(V + ( z oy) rMz ^l ( 28 ) 
Pre-multiplying (25) by (fT) 7 " , where 1 < v £ N and v * £ , and using (1), (2), and (26) yields 

-^Mz; 2)/( . - 4Mzf 2) .. - ^Mzf,, - ^Mz^ - ^ ] 
By (26), (28) and (29) we have determined z\ m . 

[73] We proceed now to derive the perturbation solutions for the general p-th 
order terms in (5) and (6). Equating the coefficients of the 5G,8Gj — SG,5G, terms of order p 
in (14) yields 
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+ 2!(p - 2)![^ 2) ,Mzf p _ 2) ,..,, + 4 )(; Mzf p . 2) ...... +- + ^Mzf^.J + - + p!^.^M#* 

Expanding z^. .„ using normalized eigenvectors of (1) as basis vectors, one has 

^ P y,.,=t^J U ( 31 ) 
where P ( k p)iJ ... stu is the coefficient of the w-th term and the number in the parenthesis in its 
subscript corresponds to that of z k (p)iJ ,.. sr Pre-multiplying (30) by {</> k ) T and using (1), (10), 
(31) and orthonormality relations of eigenvectors yields 



=^^) r t(P-0!(^<.,)>.,, + fr z U + "- + f| z U--*> 

- ( P - DK^yM^H, , ... + W^,,, + • • • + ^Mzf,.,),.., ) (32) 

- 2Kp - 2)!(^ 2K/ mz; p . 2)/ ..., + ^mz;^,.., + ... + iV Mz U,.,) — • 

-(p - 2)!2!(^. 2) ,..,Mz; 2) .. + ^. 2)> -.« Mz W" + • • • + - Mz (V<M 
The j!>-th order sensitivities of eigenvalues are obtained from (7) and (32). They depend on 
the eigenvalue and eigenvector sensitivities of orders up to p-2 and p-\ respectively. 
Substituting (31) into (13) and using (2) yields 



+ 2 \{p - 2)! [(zf 2) .. f Mz^.., + (z ( * 2) ,, ) r Mz; p . 2) .... j +--- + (z; 2>! ,) r Mz; p . 2) .,..J + - 
+ (p - 2)! 2! [(z; p . 2) ,.. j( ) r Mz[ 2) .. + (z k {p _ 2) ) r Mz; 2) , + • • • + (zf p . 2) ,... r f Mz^ ] 
+ (p - 1)! [(z[ p . l)y ,... M ) r Mz[ 1)( + (z^...,, ) r Mzf lV + ■ • • + (z;,.,).,.., f Mz; i)( ]} 



(33) 



Pre-multiplying (30) by (f) T , where 1 <v< TV and v** , and using (1), (2), and (31) yields 
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<'»~**~ p\{X l ' LKy '"dG^J-" dG.-w* dG, 

-(p-l)!(^ )f Mz^ 1)7 ,.. J( + ^ v Mz; p . t) ,.., +"' + 4MzU,,) (34) 

- 2!(/> - 2)!(^ 2) ,Mz^ 2)/ ..., ( + ^Mz^,., + • • • + ^MzJ^..., ) 

- - (P - D!(^.,)/ « Mz OV + ^W-^OV + ' ' * + 4W-. Mz 0V » 

By (31), (33) and (34) we have determined z k (p)IJ ...„. The p-th order sensitivities of 

eigenvectors can be subsequently obtained from (7). They depend on the eigenvalue and 
eigenvector sensitivities of orders up to p-\. 

[74] Equations (5) and (6) serve both the forward and inverse problems. In the 
former one determines the changes in the eigenparameters with changes in the stiffness 
parameters. Damage detection is treated as an inverse problem, in which one identifies 
iteratively the changes in the stiffness parameters from a selected set of measured 
eigenparameters of the damaged structure: X k d (k = 1,2,- -,« A ) and (k = 1,2,- •-,#!,), where 
\<n x ,n^<N . Here X k d and are assumed to be the perfect eigenparameters; simulated 
noise is included in the measured eigenparameters in the beam and frame examples that 
follow. Often we choose a set of n measured eigenparameter pairs to detect damage, i.e., 
n A =n^=n. Let the number of the measured degrees of freedom of f>/ be N m ; N m = N and 
N m < N when we have complete and incomplete eigenvector measurements, respectively. 
With reduced measurements the unmeasured degrees of freedom of {#/ is estimated first 
using a modified eigenvector expansion method (see the beam and frame examples below) 
and fa k is mass-normalized subsequently. Only the component equations corresponding to 
the measured degrees of freedom of <j> d k are used in (6). The system equations in Eqs. (5) 
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and (6) involves n x -¥n <f> N m scalar equations with m unknowns, which are in general 
determinate if n x +n^N m = m, under-determined if n x +njN m <m, and over-determined if 
n x +n <j> N m >m. In the first iteration, X k and <f> k in (5) and (6) correspond to the 
eigenparameters of the structure with the initial stiffness parameters Gj 0) . With the 
perturbation terms determined as shown above, the changes in the stiffness parameters 
SG?\ where the number in the superscript denotes the iteration number, can be solved 
from (5) and (6) using an optimization method discussed below. The estimated stiffness 
parameters of the damaged structure are updated by G, 0) = G< 0) +8Gf ) . With the updated 
stiffness parameters, the eigenparameters, X k and <j> k , in (5) and (6) are re-calculated from 
the eigenvalue problem (1) and the perturbation terms are re-evaluated. One subsequently 
finds SG\ 2) using the same method as that in the first iteration. This process is continued 
until the termination criterion, \SG\ L) \ < e , where L is the last iteration number and s is 
some small constant, is satisfied for all i= 1,2, Note that after the #>-th (\<w<L) 
iteration, we set G\ w) to G hi if G\ w) > G hi , and to zero or some small stiffness value s G if 
G\ w) <0. A flowchart for the iterative algorithm is shown in Figure IB. When a single 

iteration is used, the iterative method becomes a direct method. 

[75] Figure IB shows steps included in a method for detecting structural damage in 
accordance with one embodiment of the present invention. The method includes as an 
initial step measuring one or more eigenparameters, A*, <j) k d (Block 1). These 
eigenparameters are then compared with estimated eigenparameters associated with the 
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stiffness parameters, G, (0) (Block 2). The differences between the measured and estimated 
eigenparameters are then used by the perturbation method to establish system equations (5) 
and (6) (Block 3). The perturbation method may be a first or higher order multiple- 
parameter perturbation procedure. 

[76] Next, an optimization method may be used to find the changes in the 
stiffness parameters 8G\ w) (Block 4). These values are then compared to the predetermined 
value s (Block 5), and if the absolute values are less than s the stiffness parameters 
identified are set as G\ w ' x) (Block 6). 

[77] If the absolute values are not all less than s , the process proceeds along an 
iterative path where the stiffness parameters are first updated (Block 7). The stiffness 
parameters are then bounded between 0 or s G and G hi (Block 8), and eigenparameters 
associated with the updated stiffness parameters are calculated (Block 9). The comparison 
of Block 2 is then performed based on these calculated, or estimated, eigenparameters. 

Optimization Methods 

[78] Neglecting the residuals of order p+1 in (5) and (6) yields a system of 
polynomial equations of order p. When n A -hn^N m <m, SG\ w) (/ = 1,2, -" 9 m) at the w-th 
iteration can be solved from the resulting equations. There are generally an infinite number 
of solutions when n A + n <f) N m < m , a unique solution when n x + n (f) N m = m and p = l 9 and a 
finite number of solutions when n x +n^N m =m and p>\. When n^+n^N m >m, one 
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generally cannot find £G/ w) to satisfy all the equations, and an optimization method can be 
used to find SG\ w) which minimize an objective function related to the errors in satisfying 
these equations. We use here the same notations, e\ and e£, to denote the errors in 
satisfying the system equations in (5) and (6) respectively. Consider the objective function 

where W k x (k = l,2,---,w A ) and Wf (k = 1,2, — , w,) are the weighting factors, and / is a 
function of SGj w) when one substitutes the expressions for e\ and e£ in (5) and (6) into 
(35). When the first-order perturbations are retained in (5) and (6), the least-squares method 
can be used to determine SG\ w) which minimize (35) with W k x =Wf =1 . Other weighting 

factors can be handled by dividing (5) and (6) by / and / respectively. The method 

determines essentially the generalized inverse of the resulting system matrix, and is also 
known as the generalized inverse method. When the perturbations up to the p-th (p^\) 
order are included in (5) and (6), the gradient and quasi-Newton methods [25] can be used to 
determine SG\ w) iteratively. Unlike the generalized inverse method, the methods are 
applicable when other objective functions are defined. While the optimization methods are 
introduced for over-determined systems, they can be used when n x + n^N m < m , in which 

case J = 0 (i.e., e x = = 0) when the optimal solutions are reached. 
Gradient Method 
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[79] To minimize the objective function in (35) at the ^-th iteration, one can use 
the algorithm 

*Gg=*GW,-a»f w (36) 
to update the changes in the stiffness parameters, where ^G^ = (SG$ ) ,SG$ ) ,"',SG%j t) ) T , 

a b >0 is the step size, and f 6 _, equals the gradient vector g 6 _, = (-^,-^ r ,...,-|^.) r 

associated with SG[ b } x) . Note that the subscript b (b > 1) in all the variables in (36) denotes 
the number of nested iterations. The initial values used are SG-fy = 0 . The nested iteration 
is terminated when a b \g b _ x \\ w < y , where || • ^ is the infinity norm and y is some small 
constant, or the number of nested iterations exceeds an acceptable number, D . 



Quasi-Newton Methods 

[80] Due to its successive linear approximations to the objective function, the 
gradient method can progress slowly when approaching a stationary point. The quasi- 
Newton methods can provide a remedy to the problem by using essentially quadratic 
approximations to the objective function near the stationary point. The iteration scheme of 
these methods is given by (36) with i b _ x = B 6M g w , where B b _ x is an approximation to the 
inverse of the Hessian matrix used at the b-th nested iteration, and the other variables the 
same as those previously discussed. Initially, we set SGj$ = 0 and B 0 = I , the identity 
matrix. The matrix B b is updated using either the DFP formula 

B _ B . (^ig-^goK^ig-^a,) 1, [B t . l (g ft - gft . 1 )][B ft „(g t -g t „)r 
* (^G!; ) ) -JG£ ) ) r (g i -g 4 . l ) (g i -g 4 - 1 ) r B 6 .,(g 4 -g 6 . I ) 1 ) 
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or the BFGS formula 

= " w +U (*Gg? - *Gffl,) r (b -^) J (*G« -*Gg>/ (fc -8.-.) (3g) 

(<5Gg-^ 

The nested iteration is terminated when flfjB w g w |L<y or the number of iterations 
exceeds D. A flowchart for the quasi-Newton methods, including the step size search 
procedure as outlined below, is shown in Figure 2. 

Step Size Search Procedure 

[81] The optimal step size for the gradient and quasi-Newton methods is 
determined in each nested iteration to minimize the function J(SG\ b l l) ~-a b f b _ l ) = F(a b ) 
with respect to a b . The search procedure is divided into two phases: an initial search to 
bracket the optimum a b and a golden section search to locate a h within the bracket, as 
shown in Figure 2. 

[82] Initial Bracketing. Choose the starting point x x = 0 and an initial increment 
A > 0 . Let x 2 = jc, + A , F } = F(x x ) , and F 2 = F(x 2 ) . Since for the gradient and quasi- 
Newton methods, f 0 = g 0 and it is along a descent direction of / when A is sufficiently 
small, one has F 2 < F r Rename 2 A as A , and let x 3 = x 2 + A andF 3 = F(x 3 ) . If F 3 >F 2 , 
stop and a b is contained in the interval (*,,x 3 ). Otherwise, rename x 2 as x x and x 3 as x 2 , 
then F 2 becomes F x and F 3 becomes F 2 . Rename 2 A as A , and let x 3 = x 2 + A and 
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F 3 = F(Xj). Compare F 3 andF 2 , and repeat the above procedure if F 3 < F 2 untilF, >F 2 , 
with the final interval (x,,x 3 ) containing cc b . 

Golden Section Search. If |x 3 - x 2 \ > \x 2 - x, \ , define a new point: 

x 4 =x 2 +0.382(x 3 -x 2 ) (39) 
Otherwise, rename x, as x 3 and x 3 as jc„ and then define x 4 using (40). Let F 4 = F(x 4 ). If 
F 2 <F 4 , rename x 4 as x 3 , then F 4 becomes F 3 . Otherwise, rename x 2 as x, and x 4 as x 2 , 
then F 2 becomes F, and F 4 becomes F 2 . Compare |x 3 -x 2 | and |x 2 -x,|, and repeat the 
above procedure until \x 3 -x x \ |f w || w < s. , where e a is some small constant. Then choose 



a. =— - 

* 2 



Mass-Spring System Example 

[83] The algorithm discussed above is used to identify the stiffness parameters of a 
N-degree-of-freedom system consisting of a serial chain of masses and springs, such as the 
system shown in Figure 3. Let the masses of the system be M ( . =1 kg (i = 1,2, — ,N), and 
the stiffnesses of the undamaged springs be G hi =\Nlm (i = 1, 2, —,m), where m = N + \. 
The system is said to have a small, medium and large level of damage if the maximum 
reduction in the stiffnesses is within 30%, between 30 and 70%, and over 70%, respectively. 
The mass matrix M is an JVxJV identity matrix, and the stiffness matrix K is a banded 
matrix with entries K,. = G i + G M (i - 1,2,-, N), K f(f+1) = K {MV = ~G M (i = 1,2,-, N -1), 

and all other entries equal to zero. The matrices —and — - have a unit value in entries 

OLr x OLr N 
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(1, 1) and (N, N), respectively, and zero entries elsewhere. The nonzero entries of the 
matoces — (, = 2,3,-,tf-l)are (—) ( m X m, =(—)*= 1 and (— ), (W) = (— ) (MV = -1. 



dG i 



dG : 



dG/ 



8G ; 



[84] We look at a forward problem first with N = 3 and m = 4 . The stiffnesses of 
the damaged system are = G rf3 = 1 Nlm , G d2 =03N/m and G d4 =0N/m. The 
undamaged system is considered as the unperturbed system and the damaged system as the 
perturbed system. Based on the eigenparameters of the undamaged system, the 
eigensolutions of the damaged system are obtained using the first-, second-, and third-order 
perturbations, as shown in Table 1 below. The results show that even with the large changes 
in stiffness, the third-order perturbation solutions compare favorably with the exact 
solutions for the damaged system. The higher-order perturbation solutions can be used for 
large order systems when their direct eigensolutions become costly. 

Table 1. Eigensolutions of the damaged system from an eigenvalue problem solver 
(exact) and perturbation analysis 



Eigenparameters 



Exact 



Unperturbed 



Perturbed 



1st order 



2nd order 



3rd order 



A 2 
A 2 



0.10602 

1.27538 

2.21859 

0.16516 
0.65733 
0.73528 

0.95541 
0.07840 
-0.28470 

0.24478' 
-0.74952 
0.61507 



0.58579 

2.00000 

3.41421 

0.50000 
0.70711 
0.50000 

0.707 11 
-0.00000 
-0.70711 

0.50000 
-0.70711 
0.50000 



0.30576 

1.15000 

2.14424 

0.28523' 
0.68836 
0.74129 

0.95459' 
0.10607 
-0.45962 

0.36477 
-0.72586 
0.60871 



0.15670 

1.25500 

2.18830 

0.161 U 
0.66444 
0.79452 

' 1.00188] 
0. 1 43 1 9 [ 
-031776} 

0.29635 
-0.74132 
0.62482 



0.10090 

1.29344 

2.20567 

0.12907' 
0.65040 
0.76653 

0.97976' 
0.12310 
-0.26810 

0.26445 
-0.74875 
0.62362 
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[85] Consider now the damage detection problem with N = 9 y /w = 10, 
G d5 =0.5N/m y G dz =0JN/m, G dl0 =0.8N/m and G di =lN/m (/ = 1,2,3,4,6,7,9). We 
set W* = Wf = 1 , s = 0.001 , r = 10" 10 , a { = 1 for all /, and Z) = 500; the actual numbers of 
nested iterations in all the cases are much smaller than D. Since vanishing stiffness in any 
spring other than the two end ones in Figure 1 can result in two decoupled subsystems, 
whenG^ <0 we set G\ w) to £ G =0.\N/m in the first two iterations and to 0.01 N/m in the 
remaining iterations for all the cases considered here. A relatively large value is assigned to 
s G in the initial iterations to avoid close eigenvalues in the mass-spring system and small 
denominators in (19). This improves convergence especially when a small number of 
eigenparameter pairs are used. A smaller value is used for s G in the later iterations to 
improve the accuracy of stiffness estimation when there is a large level of stiffness reduction. 
Using the first-order perturbations and different numbers of eigenparameter pairs, the 
maximum errors in estimating the stiffnesses of the damaged system at the #>-th iteration, 
defined by 

E = max J L (40) 

are shown in Figures 4A and 4B for all the iterations. In Figure 4A p = 1 and n = 1,2,3, and 
in Figure 4B p = 1 and n = 4,5,-*-,9. When « = 1, the error decreases slowly, though 
monotonically, and there is an estimation error of 1.5% at the end of iteration. While the 
errors can increase with the iteration number before approaching zero for «=3, they 
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decrease monotonically for n = 2 and n > 4 . All the stiffnesses are exacdy identified at the 
end of iteration when n > 2 . Note that the number of the system equations equals and 
exceeds the number of unknowns when n = 1 and n > 2 , respectively. Since the system 
equations are linear, they have a unique solution when n = 1 , and / has a unique minimum 
when n > 2 . With the small y the gradient method and the quasi-Newton methods using 
the DFP and BFGS formulas yield exactly the same results as the generalized inverse 
method (not shown here). Because the generalized inverse method does not involve any 
nested iteration, it is the most efficient one among the four methods. While not shown here, 
the results indicated that that the quasi-Newton methods converge faster than the gradient 
method and the BFGS method has the similar performance to the DFP method. In what 
follows the BFGS method will be used with the higher-order perturbations. With the 
second-order perturbations the errors shown in Figure 4C decrease monotonically for all n 
(n = 1,2, ••■,9). The errors atw = l in the expanded view decrease in the order 
n = 3,2,4,9,7,8 , with the lines for n = 5 and 7 virtually indistinguishable. In Figure 4(A-C) 
the following symbols are used for different n\ * « = 1; -e— n = 2; A n = 3; 
-*~„ = 4', n = 5; « = 6; n = 7 ; * n = S; ° * = 9 . While the use 

of the second-order perturbations improves the accuracy of stiffness estimation in each 
iteration and reduces the number of iterations, it takes a much longer time to compute the 
higher-order perturbations and the associated optimal solutions. 

[86] When only the first few eigenvalues are used, for instance, n x = 5 and n 4 = 0 , 

the stiffnesses identified with the first-order perturbations, 



32 



G (4) = (0.875,0.976,0.926,0.864,0.699,0.699,0.864, 0.926, 0.976, 0.875) r Nlm (41) 

where the number in the superscript denotes the last iteration number, correspond to those 
of a different system with the same eigenvalues for the first five modes as the damaged 
system. The same stiffnesses are identified with the second-order perturbations. Similarly, 
when the first eigenvector is used, i.e., n $ = 1 and n k = 0, the stiffnesses identified with the 

first-order perturbations are those of a different system with the same eigenvector for the 
first mode as the damaged system: 

G (9) = (0.989, 0.991, 0.995, 1, 0.520, 0.829, 0.940, 0.668, 0.959, 0.768) r Nlm (42) 
With the second-order perturbations the stiffnesses of the damaged system are identified. 
The stiffnesses identified are not unique because the system equations in each iteration are 
under-determined. The solution given by the generalized inverse method here in each 
iteration is the minimum norm solution. Increasing the number of eigenparameters used can 
avoid this problem. 

[87] If the system has a large level of damage, i.e., G ds = 0.3 jV I m , G dx0 = 0.1 N I m , 

with the other parameters unchanged, the stiffnesses of the damaged system are identified 
with the first-order perturbations after 55 iterations when n = 1 and 6 iterations when n = 2 
and 3 as shown in Figure 5A. For n = 4,5,6,7, the errors decrease monotonically and the 
number of iterations is reduced slightly, as shown in Figure 5B, which has an expanded view 
near w = 1 . With the second-order perturbations the errors shown in Figure 5C decrease 
monotonically for n = l,2,---,5 5 and the number of iteration for n = 1 is reduced from 55 in 
Figure 5A to 4. The errors at w = 1 in the expanded view in Figure 5C decrease in the order 



33 



n = 3,2,4,5 , with the lines for n = 2 and 4 virtually indistinguishable. The symbols used in 
Figure 5A-C for n = 1, 2, • • • , 7 are the same as those in Figure 4A-C 

[88] Finally, consider a large order system with a large level of damage: N = 39 , 
m = 40, G rfl2 =0.77V/m, G, I9 = G rf37 = O.ltf/m , G, 28 = 0.8A^/m , G di =lN/m (l<z<40 
and 12,19,28,37), and the other parameters are the same as those in the previous 
example. With the first-order perturbations the exact stiffnesses are identified after 57 
iterations when n = 1 , as shown Figure 6A. Using a larger number of eigenparameter pairs 
(« = 2,3,4,5) significantly reduces the number of iterations, as shown in Figure 6B. The 
errors at w = 2 in the expanded view in Figure 6B decrease in the order n = 4,2,3,5 . With 
the second-order perturbations the exact stiffnesses are identified after 6 iterations when 
n = 1 and 3 iterations when n = 2 , as shown in Figure 6C, which has an expanded view for 
l<w<4. The symbols used in Figure 6(A-C) for « = l,2, — ,5 are the same as those in 
Figures 4A-C and 5A-C. 

Fixed-Fixed Beam Example 

[89] The algorithm discussed above may be applied to detecting structural damage 
in an aluminum beam with fixed boundaries. The beam of length L t =0,7/n, width 
W = 0.0254 m, and thickness H = 0.0031 m has an area moment of inertia 

/ = -^WH 2 ' = 6.3058 x 10" n m A and a mass density p = 2715 kg/m 3 . The finite element method 

is used to model its transverse vibration. The beam is divided into N e elements, as shown in 
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Figure 7, with the length of each element being /. = ^. There are iV e +1 nodes. With V, 
and B, denoting the translational and rotational displacements at node i (i = l,2,-,# e +1), 
the displacement vector of the i -th (i = 1,2,-,AT.) element is M,F, +1 ,0, +1 f . The Young's 
modulus is assumed to be constant over each beam element and that of the /-th element is 
denoted by G, . The Young's modulus of the undamaged beam is G h =69xl0 9 N I m 1 . Hence 
G =G for i = l,2,-,m, where ro = A/ e . Small to large levels of damage correspond to 

hi h 

reductions of the moduli defined for the Mass-Spring Example. The mass and stiffness 
matrices of the / -th beam element are 



420 



156 

-22/, 

54 



-22/ e 
-13/ 



54 13/, ] 

-13/, -3/ e 2 
156 22/ 



-3/ 2 22/ 4/ 2 



Kf = 



GJ 



12 

6/. 
-12 

61. 



61, -12 

All -61. 

-6l e 12 

111 -61. 



6/ e 
211 

-6/ e 
Ml 



(43) 



Using the standard assembly process yields the 2(N e + 1) x 2(N e + 1) global mass and stiffness 
matrices. Constraining the translational and rotational displacements of the two nodes at the 
boundaries to zero yields the NxN M and K matrices, where N = 2(N e -1) is the degrees 
of freedom of the system. The displacement vector of the system, involving the 
displacements of the 2 nd through N e -th node, is [K 2 ,^ 2 ,F„5„-,K W> ,^J . The matrix — 

(i = l ( 2,--,m) can be obtained from K by setting G i =1 and G, = - =G M =G M =-- = G N =0. 
The parameters W* , W} , s , y , <T <} andD are set to the same values as those previously 
discussed, and e c is set to 0. 1 5G h in the first two iterations and to0.05G A in the remaining 
iterations. The first-order perturbations are used unless indicated otherwise. 
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[90] Consider first the cases with N e =m = \0 and N = 18 . When the system has a 
medium level of damage: 

G, = (1, 1, 1, 1, 0.5, 1,1, 0.7, 1, 0.8) r x G h (44) 
the stiffness parameters of the damaged system are identified after 6 iterations with n - 1 . 
When the system has a large level of damage: 

G d = (1, 1, 1, 1, 0.3, 1, 1, 0.7, 1, 0. If x G h (45) 
the stiffness parameters of the damaged system are identified after 7 iterations with n = 1 . 
Consider next the cases with N e = 20, 40 and 80. For the systems with medium and large 
levels of damage, the stiffness parameters of the first 10 elements are given by (44) and (45), 
respectively, and those of the remaining elements are G h . In all the cases the stiffness 
parameters of the damaged systems are identified within 10 iterations when n = \. The 
numbers of iterations are reduced slighdy when the second-order perturbations are used. 
Note that the system equations are over-determined when n = 1 . 

[91] When only the translational degrees of freedom of an eigenvector are 
measured, a modified eigenvector expansion method is used to estimate the unmeasured 
rotational degrees of freedom. To this end, </> d is partitioned in the form 

$ = [(^<L) f > (0duYf > where <j) k dm and <j> k du are the measured and unmeasured degrees of 
freedom of <f> k dy respectively. Similarly, </> k in (6) is partitioned in the form 
<l> k =[(<fit) T > (fiuYf , where <j> k m and <f> k correspond to the measured and unmeasured 
components of <f> d , respectively. Since <t> dm , <f) k m and <j) k are known in each iteration, <f> du is 
estimated from (/> du =[(#) + ^J#, where the superscript + denotes generalized inverse. 
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Once the rotational degrees of freedom of <f) k d are determined, (j> k d and ^*are converted to 
their original forms and <j> k d is mass-normalized. Only the component equations 
corresponding to the measured degrees of freedom of $ are used in (6) and the system 
equations are determinate when n = \. The exact stiffness parameters of the damaged 
systems considered above can be identified. For the 10- and 20-element beams with the 
medium levels of damage, the stiffness parameters are identified after 6 and 24 iterations, 
respectively, when n x = 1 and n 4 = 2 . For the 10- and 20-element beams with the large levels 

of damage, the stiffness parameters are identified after 9 iterations when n A - 1 and n 4 - 3 
and 10 iterations when n = 2 , respectively. 

[92] Finally, the effects of measurement noise on the performance of the algorithm 
are evaluated for the 10-element beam with the large level of damage. Simulated noise is 
included in the measured eigenparameters: 

A* = K k + vR'X . = € + o*lff (46) 

where A*/ and f k are the k -th perfect eigenvalue and eigenvector, respectively, R k is a 
uniformly distributed random variable in the interval [-1,1], is a diagonal matrix whose 
diagonal entries are independently, uniformly distributed random variables in the interval 
[-1,1], and [0,1] is the noise level. Note that R k x and R* are generated for each 
measured mode. Each random parameter is generated 10 times and the average is used. 
Three different noise levels are considered: u = 5%, 10% and 20%. When all the degrees of 
freedom of an eigenvector are measured, the stiffness parameters identified with n = 1 , 2, 
and 3 are shown in Figure 8A, 8B, and 8C, respectively. The following symbols are used for 
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different noise levels: S , o = 5% ; B , u = 10% ; B , o = 20% . When only the translational 
degrees of freedom of an eigenvector are measured, the eigenvector expansion method 
described above is used and the stiffness parameters identified with n = 2 and n = 3 are 
shown in Figure 9A and 9B, respectively. The same symbols as those in Figure 8A-C are 
used in Figure 9A-B for different noise levels. The stiffness parameters corresponding to 
u = 0 in Figures 8 and 9 are the exact values. It is seen that in the presence of noise, the 
stiffness parameters can be accurately identified with an increased number of measured 
eigenparameters. 



SPACE FRAME EXAMPLE 

[93] The damage detection method can be applied to more complex structures, 
such as, for example, the modular, four bay space frame shown in Figure 10. In this 
example and all the following examples the first order perturbation approach along with the 
generalized inverse method is used. The four-bay space frame example shown in Fig. 10 is 
made of extruded aluminum with 48 members, and is 8' tall, 1.46' wide and 1.8' deep. The 
horizontal and diagonal members have the same cross-sectional dimensions 
(25.4 mm x 12.7mm) and the vertical members are "L"-angles with cross-sectional dimensions 
50.8mm x 50.8mm x 6.35mm (thickness). All the members are connected through bolted joints. 
The frame is assumed to be fixed to the ground. The finite element method is used to 
model the 3-dimensional vibration of the frame, with each member modeled with four 12- 
degrees-of-freedom beam elements. The total degrees of freedom are 960 when the 
boundary conditions are applied. The Young's modulus is assumed to be constant over each 
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member and that of the z-th (1</</w=48) member is denoted by G r The Young's 
modulus of the undamaged member is G h =69xl0 9 N I m 1 and G hi =G h . A vertical (member 
1) and a diagonal (member 48) member in the first bay are assumed to have 70% and 30% 
reductions in Young's modulus, respectively, and all the other members are assumed to be 
undamaged. The first two vibration modes, corresponding to the bending of the frame 
along the x and y directions, respectively, are used to detect damage. The translational 
degrees of freedom of the 16 nodes, denoted by "S" in Figure 10, along the x and y 
directions are assumed to be measured with 10% measurement noise. All the other degrees 
of freedom are estimated using the eigenvector expansion method discussed above and the 
measured modes are mass-normalized. The Young's moduli of all the members of the 
damaged truss are identified as shown in Figure 10, with the maximum estimation error less 
than 7%. Note that the truss has closely spaced vibration modes. With the modes arranged 
in the order of increasing frequencies in each iteration, the method has effectively handled 
mode switching that has occurred in the damage detection process. 

[94] In summary, the damage detection method identifies stiffness parameters in 
structures, which have a small, medium, and large level of damage if the maximum 
reduction in the stiffnesses is within 30%, between 30 and 70%, and over 70%, respectively. 
A large level of damage is studied in many examples because this poses the most challenging 
case, with sever mismatch between the eigenparameters of the damaged and undamaged 
structures. 

[95] The damage detection method as embodied and broadly described herein can 
be applied to structures that can be modeled with beam elements. A beam element is an 
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element that has one dimension that is much longer than the other two. This element is very 
good at modeling 'T'-beams, rectangular beams, circular beams, "L"-angles, "C"-channels, 
pipes, and beams with varying cross sections. Structures that can be modeled with this 
element include, but are not limited to, lightning masts, light poles, traffic control poles, 
pillar type supports, bridges, pipelines, steel building frameworks, television, radio, and 
cellular towers, space structures, cranes, pipelines, railway tracks, and vehicle frames. 
Structures that can be modeled as beam elements are used simply for ease of discussion, and 
it is well understood that the damage detection method discussed above may be applied to 
structures that can be modeled by other elements using the finite element method or 
modeled by using other methods. 

Damage Detection Using Changes of Natural Frequencies: Simulation and 
Experimental Validation 

[96] For structures such as beams and lightning masts in electric substations, using 
only the changes in the natural frequencies can relatively accurately detect the location(s) and 
extent of damage, even though the system equations are severely underdetermined in each 
iteration. This is an interesting finding as it is much easier to measure the natural frequencies 
than the mode shapes, and demonstrates the effectiveness of the iterative algorithm. 
Extensive numerical simulations on beams and lightning masts confirmed this finding. 
Experiments on the beam test specimens with different damage scenarios and a lightning 
mast in an electric substation validated the simulation results. The beam test specimens and 
the lightning mast are used as examples for demonstration purposes, and the method can be 
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applied to other structures. Note that unlike the beam example shown earlier, where the 
Euler-Bernoulli beam finite element model is used, the Timoshenko beam finite element 
model is used in all the examples here. The Timoshenko beam theory is found to be more 
accurate in predicting the natural frequencies of the lightning masts and circular beams than 
the Euler-Bernoulli beam theory. 

[97] For a cantilever beam, simulation results show that the damage located at a 
position within 0-35% and 50-95% of the length of the beam from the cantilevered end can 
be easily detected with less than 5 measured natural frequencies, and the damage located at a 
position within 35-50% of the length of the beam from the cantilevered end and at a 
position within 5% of the length of the beam from the free end can be relatively accurately 
detected with 10-15 measured natural frequencies. 

Numerical and Experimental Verification 

Cantilever aluminum beams 

[98] Experimental damage detection results for four different scenarios are shown 
first, followed by various simulation results. 

❖ Scenario 1: Evenly-distributed damage machined from the top and the 
bottom surfaces of the beam test specimen. 

[99] The aluminum beam test specimen shown in Figure 12 is 45 cm long by 2.54 
cm wide by 0.635 cm thick. It is divided into 40 elements (each element has a length of 
1.125 cm). The beam has a section (from approximately 10 cm to 15 cm from the 
cantilevered end) of 5 cm long and 7.62E-4 m thick machined both from the top and the 
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bottom surfaces of the beam. This corresponds to 56% of damage (or reduction of bending 
stiffness El) along the length of five elements (from the 9* to the 13 rd element). Using the 
changes of the first 2 to 5 measured natural frequencies, damage is detected within 7 
elements using 2 or 5 measured frequencies (from the 7* to the 13 rd element with 5 
measured frequencies and from the 5 th to the 11* element with 2 measured frequencies) and 
within 8 elements using 3 or 4 measured natural frequencies (both from the 6 th to the 13 rd 
element), as shown in Figure 13. With 5 measured frequencies the average damage of the 7 
elements in Figure 13, from 6.75 cm to 14.625 cm, is 46%. Note that the elements with 
damage less than 10% are not accounted for here. The extent of damage detected is slightly 
lower than the actual extent because the predicted damage occurs at 2 more elements (the 7 th 
and 8 th elements) than the actual one. The error results from the solution of the severely 
underdetermined system equations (5 equations with 80 unknowns). 

❖ Scenario 2: The same aluminum beam test specimen as in Scenario 1 
clamped at the other end. 

[100] The same beam was tested in the clamped- free configuration with the clamped 
end reversed, which placed the damage from 25 cm to 30 cm (from the 23 rd to the 27 th 
element) from the cantilevered end. The damage detection results with 3 to 5 measured 
frequencies are shown in Figure 14. With 5 measured frequencies the average damage 
between 23.625 cm to 32.625 cm (from the 22 nd to the 29 th element) is 40%. Again the 
elements with damage less than 10% are not included. 

❖ Scenario 3: Undamaged cantilever aluminum beam test specimen with the. 
same dimensions as above. 
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[101] An undamaged aluminum beam test specimen was clamped at one end with 
the same configuration as shown in Figure 12. With the first 3 to 4 measured frequencies 
the maximum error for the estimated bending stiffnesses of all the elements is within 10%, 
as shown in Figure 15. 

❖ Scenario 4: A cut of small width on a cantilever aluminum beam test 
specimen, shown in Figure 16 with the same dimensions as above. 

[102] The beam shown below has a cut that is 0.4191 cm deep and .1016 cm wide, 
which corresponds to a 96% reduction in bending stiffness at the cut. The beam is divided 
into 45 elements and the cut is located in the middle of the 23 rd element. With 3 to 5 
measured natural frequencies, the damage is detected at several elements surrounding the 
23 rd element, as shown in Figure 17. It was found there was a roughly 50%, 40%, and 60% 
reduction in stiffness at the 22 nd , 23 rd , and 24 th element, respectively. The estimated damage 
extent at the above elements is lower than that at the cut because the damage is distributed 
along these elements. 

[103] To examine the effectiveness and robustness of the damage detection 
algorithm, various simulations with different damage scenarios were carried out. In this way, 
we can gain more insight concerning the accuracy of the finite element model, convergence 
of the estimated bending stiffnesses of all the elements of the beam with the increased 
numbers of measured natural frequencies and/or mode shapes, and region of the beam 
within which the damage can be detected with few measured frequencies. With the 
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cantilever aluminum beam divided into 40 elements, the following two simulations have the 
similar damage location and extent to those in Scenarios 1 and 2 in the experiments: 

❖ Simulation 1: Uniform damage between 9.0 cm and 15.75 cm from the 
cantilevered end of the beam with the same dimensions as discussed above. 

❖ Simulation 2: Uniform damage between 29.25 cm and 36 cm from the 
cantilevered end of the beam with the same dimensions as discussed above. 

[104] Simulation results in Figure 18 and Figure 19 show that the damage can be 
relatively accurately detected with the first 3 to 5 measured frequencies, and can be 
accurately detected with the first 10 measured frequencies. 

[105] With the cantilever aluminum beam divided into the same number of 
elements, two more simulations are presented here: one for a multiple damage scenario - 
70% of damage at the 3 rd element and 30% of damage at the 20 th element, and the other for 
a 50% of uniform damage from the 16 th to the 18 th element (i.e., the damage is located at a 
position within 35-50% of the length of the beam from the cantilevered end). 

❖ Simulation 3: Multiple damage of the beam with the same dimensions as 
discussed above: 70% of damage at the 3 rd element and 30% of damage at the 
20 th element. 

[106] Simulation results in Figure 20 show that the multiple damage can be relatively 
accurately detected with the first 3 to 5 measured frequencies, and accurately detected with 
the first 10 measured frequencies. 

❖ Simulation 4: Uniform damage from the 16 th to the 18 th element of the beam 
with the same dimensions as discussed above. 
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[107] Simulation results in Figure 21 show that the damage within 35-50% of the 
length of the beam from the cantilever end can be accurately detected with 10-15 measured 
frequencies. 

Lightning Masts 

[108] The lighting mast shown in Figure 22 has two sections of constant cross- 
sections and an eccentric spike. The lengths of the lower and upper sections are 6.89 m and 
6.83 m, respectively, and that of the spike above the upper section is 1.4375 m. Both mode 
shapes and natural frequencies were measured for this mast. The mode shapes were 
measured by using a laser Doppler vibrometer. The natural frequencies were measured 
using both the laser Doppler vibrometer and an accelerometer. It takes about .5-1.5 days to 
measure the mode shapes and about 30 minutes to measure the natural frequencies. It 
would be even more difficult to measure the mode shapes for some of the taller lightning 
masts, which are 100 and 130 feet tall. A finite element model was made using OpenFEM. 
Once the model was completed the measured and calculated natural frequencies and mode 
shapes were compared. The mast was expected to be undamaged, since by inspection the 
measured and calculated natural frequencies matched (Table 2). It was found that the first 
mode shape measurement is affected by wind and high modes are less affected by the wind. 
The measured frequencies are not affected by the wind and can be used for damage 
detection. 

[109] Table 2. Comparison of measured and calculated (from the finite element 
model) natural frequencies 
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Mode# jMeasiredl 


FEM I Error % 


1 


1.17 


1.18 


-0.8648 


2 


5.10 


5.23 1 


-2.4104 


3 
4 


7.83 
8.17 


7.89 
8.67 


-0.8418 
-6.1104 


5 


16.60 


16.55 


0.31327 


6 


29.26 


30.36 


-3.7385 


7 


45.32 


47.85 


-5.5831 


8 


47.75 


50.87 


-6.5191 


9 


54.10 


55.81 


-3.1439 


10 


67.45 


68.83 


^.0476 


11 


74.70 


78.99 


-5.7456 



[110] Damage detection was then performed using only the first 4 to 7 measured 
natural frequencies, as shown in Table 2. The mast is modeled by 40 elements; the lower 
section has 18 elements, the upper section has 15 elements, and the spike has 7 elements. 
The experimental damage detection results are shown in Figure 23. The mast is modeled by 
40 elements with 20 elements in each of the two sections. It appears that the elements 
around the joints near the middle of the mast (8.2 m in Figure 23) and at the free end of the 
mast (13.72 m in Figure 23) had some damage. This is most likely due to the fact that the 
joints were modeled in the finite element model as being infinitely stiff, which is almost 
never the case. Also, the spike is actually connected to the free end of the mast at two 
points. In the finite element model here it is assumed that the spike is attached to the mast 
along the whole line of contact that has a length of 0.34 m, which makes the model a little 
stiffer. Similar damage detection results were obtained with different numbers of measured 
frequencies. 

[Ill] Simulations for the same lightning mast as shown in Figure 22 with different 
damage scenarios are carried out. Shown in Figure 24 are the simulation results for the mast 
with 70% damage between 0.76 m and 1.15 m from the ground. Shown in Figure 25 are the 
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simulation results for the mast with 50% damage between 6.89m and 7.35m from the 
ground. In both cases the location and extent of damage can be relatively accurately 
detected using the first 5 measured natural frequencies and accurately detected using the first 
8 or 10 measured frequencies. 

Methods to Handle Some Ill-Conditioned System Equations 

[112] When the first order perturbation approch is used, one needs to solve in each 
iteration a system of linear algebraic equations 

ASG = F (47) 

which are the linearized equations of (5) and (6), where A is the system matrix, F is the 
vector representing the differences between the measured and estimated natural frequencies 
and/or mode shapes, SG is the optimal changes in the stiffness parameters to be found. 
While the gradient and quasi-Newton methods can be used, the generalized inverse method 
is most efficient because it does not involve nested iterations, and thus SG = A + F , where A + 
is the generalized inverse of A. The problem (47) may be ill posed for certain cases, 
especially in the first several iterations, because A + can be relatively large and small changes 
in F can result in large changes in the solution SG . Since the stiffness parameters cannot be 
negative or greater than the corresponding values of the undamaged structure, the solution 
may not converge. 

[113] Ill-conditioning problems do not occur in all the examples described above. 
Sometimes they can occur. Consider, for example, a cantilever aluminum beam of the same 



47 



dimensions as those of the beam shown in Figure 12. The beam is divided into 36 elements 
and assumed to have 30% of damage at the 5 th element and 90% of damage at the 18 th 
element. The translational degrees of freedom of the first or second mode shape vector are 
used to detect damage. The system equations are determinate, and the ill-conditioning 
problem occurs in the iterations. When a natural frequency is also used in additional to an 
incomplete eigenvector described above, the ill-conditioning problem can also occur. The 
following regularization methods can be used to handle the ill-conditioning problem: 
Method 1. Estimate A + from (A T A + n*)' 1 * T , where r] is a small positive constant that can 
be searched in several ways. One way is set 7 = ^xl.618xmin(10,||^G|| m ), where denotes 
the infinity norm, with the initial value rj = \a t , until \\SG\\ a is less than 0.8, for example. 
Method 2. To constrain the magnitude of SG , we include it in the objective function to be 
minimized. For example, we can minimize the following objective function 

f(SG) = XI^- - VG,.) 2 +S(*G,) 2 (48) 

instead of (35), where N 0 =n x + n^N m is the number of equations in (5) and (6). The 
optimal solution can be obtained by using the generalized inverse method for the expanded 
system A* =[A;I] and the expanded vector F* =[F;0], where I is the wxw identity matrix 

and 0 is the mxl zero vector. 

[114] Note that the solutions from the regularization methods are not stricdy the 
optimal solutions for the original objective function in (35). With accurate and sufficient 
measurement information and proper handling of the ill-conditioning problem, the system 
equations can become well conditioned in the last few iterations and regulation does not 
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need to be applied. Consequendy the stiffness parameters can be more accurately 
determined. Sometimes regulation may over constrain the magnitude of SG , and the 
termination criterion jSGj^ <s may be satisfied during the regulation process. In this case, 

we may set £G, = ' to amplify the magnitude of SG , and the iteration is terminated 
\\ SG \l 

after the system equations become well conditioned. Since the second regulation method 
does not need to search for 77 , it can be more effective when it works. When one carries out 
the damage detection procedures using different combinations of measured frequencies and 
mode shapes and obtains the same or similar stiffness parameters, one can have sufficient 
confidence on the results obtained. 

[115] Using only the translational degrees of freedom of the first eigenvector and 
the first regulation method, the estimated bending stiffnesses of all the elements of the beam 
are shown in Figure 26A. Using the translational degrees of freedom of the second 
eigenvector and the first regulation method, the estimated bending stiffnesses of all the 
elements of the beam are shown in Figure 26B. In both cases the locations of damage are 
exactly detected and the extent is relatively accurately determined. Using the translational 
degrees of freedom of the first eigenvector along with the first natural frequency and the 
first regulation method, the bending stiffnesses of all the elements of the beam are exactly 
determined, as shown in Figure 26C. 

[116] Similarly, using the translational degrees of freedom of the first eigenvector 
and the second regulation method, the estimated bending stiffnesses of all the elements of 
the beam are shown in Figure 27A. Using the translational degrees of freedom of the 
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second eigenvector and the second regulation method, the estimated bending stiffnesses of 
all the elements of the beam are shown in Figure 27B. In both cases the locations of damage 
are exactly detected and the extent is relatively accurately determined. Using the translational 
degrees of freedom of the first eigenvector along with the first natural frequency and the 
second regulation method, the bending stiffnesses of all the elements of the beam are exactly 
determined, as shown in Figure 27 C. 

Conclusions 

[117] Thus, the sensitivities of eigenparameters of all orders may, for the first time, 
be derived using a multiple-parameter, general-order perturbation method. The higher-order 
solutions may be used to estimate the changes in the eigenparameters with large changes in 
the stiffness parameters. The perturbation method may be combined with an optimization 
method to form a robust iterative damage detection algorithm. The gradient and quasi- 
Newton methods can be used for the first or higher order system equations, and the 
generalized inverse method can be used efficiently with the first order system equations 
because it does not involve nested iterations. Including the higher-order perturbations can 
significantly reduce the number of iterations when there is a large level of damage. A 
modified eigenvector expansion method is used to estimate the unmeasured component of 
the measured mode shape. For many cases, the location(s) and extent of damage can be 
relatively accurately detected using only measured natural frequencies. Methods to handle 
ill-conditioned system equations that may occasionally arise are developed and shown to be 
effective. Numerical simulations on different structures including spring-mass systems, 
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beams, lightning masts, and frames show that with a small number of measured 
eigenparameters, the stiffness parameters of the damaged system may be accurately 
identified in all the cases considered. Experiments on the different beam test specimens and 
the lightning mast in an electric substation validated the theoretical predictions. The 
methodology can be readily applied to various operation structures of different sizes by 
incorporating their finite element models or other mathematical models. 

[118] One example of a practical application of this damage detection method is 
shown in Figure 28, in which a structure 10 may represent any type of structure that can be 
modeled by beam elements, on which periodic damage assessment must be conducted. This 
includes, but is not limited to, structures such as lightning masts, utility poles, cell towers, 
and other such structures. A structure that can be modeled by beam elements is used simply 
for ease of discussion, and it is well understood that the general order perturbation method 
discussed above may be applied to a number of different structural members without 
departing from the spirit of the invention as embodied and broadly described herein. 

[119] The structure 10 divided into elements Ei through E n , and with nodes Ni 
through N n , is equipped with a sensor 40, such as, for example, an accelerometer. The 
sensor 40 is configured to measure a response to an impact F or a force from a shaker 
applied to the structure 10. The impact F may be applied in the form of a single impact, as 
shown in Figure 28, or it may be in the form of a series of random impacts Ft through F n , as 
shown in Figure 29. The location on the beam where the impact F is applied may be varied, 
but the impact F is preferably not applied at one of the nodal points of vibration modes 
whose natural frequencies and mode shapes need to be measured. In either instance, the 
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sensor 40 senses a response in the structure 10 to the impact F or the force from the shaker, 
and transmits that response to a processor 20 equipped with, among other elements, the 
damage detection method 30 discussed above. 

[120] In accordance with the method as described above, upon receipt of the 
response signal from the sensor 40, the processor 20 accesses and performs the method 30 
as previously discussed, and, as the data converges, determines an extent and location, by 
element, of any structural damage present in the structure 10. 

[121] Although the convergence of the system to resolution is not dependent on 
where the impact F is applied to the structure 10, a signal to noise ratio of the response may 
be increased, and an efficiency of the system optimized, as the application point of the 
impact F is moved away from a fixed point No, if one needs to measure the natural 
frequencies and/or mode shapes of lower modes. Similarly, although the system will 
converge regardless of the average energy and timing of the impact F or series of impacts Fi 
- F n applied to the structure 10, a signal to noise ratio of the response may be increased, and 
an efficiency of the system optimized based on a random series of impacts or the random 
shaker excitation to the structure 10, if the structure is large or slightly nonlinear or if there is 
an ambient excitation to the structure such as wind. 

[122] The structure 10 may be excited in a number of different manners. For 
example, an impact F may be applied manually or a shaker may be used, at any given point 
on the structure 10, excluding the fixed point No or the nodal points of vibration modes 
whose natural frequencies and/or mode shapes need to be measured. However, in an effort 
to increase the signal to noise ratio in the signal provided by the sensor 40 to the processor 
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30 if the structure is large or slightly nonlinear or if there is an ambient excitation to the 
structure such as wind and to provide the convenience and portability over the shaker test, a 
series of random impacts Fi - F n shown in Figure 29 may be applied manually or by a 
specially designed device to generate the impact F. Although accurate and reliable damage 
assessments can be achieved regardless of how the impact F is applied to the structure 10, 
results may be improved using a random series impact method, which involves using a series 
of impacts F of random amplitudes and random arrival times. A series of random impacts 
has been shown to increase an energy input to the structure 10, improve the signal to noise 
ratio, especially in such situations as strong wind excitation, and average out slight 
nonlinearities that arise, for example, from bolted joints and extract linearized 
eigenparameters. 

[123] Figure 29 illustrates a structure 10, sensor 40 and processor 20 similar to those 
shown in Figure 28. An impact F in Figure 29 is applied to the structure 10, at a location 
other than No or one of the nodal points of vibration modes whose natural frequencies 
and/or mode shapes need to be measured, as a series of random impacts Fi through F n 
delivered at random amplitude and arrival time. These random impacts Fi - F n may be 
delivered manually, or, as shown in Figure 30, may be delivered in an automated fashion 
through the use of a random impact hammer device 1550. 

[124] Different methods have been employed in conventional vibration testing in 
order to excite a test specimen. Shaker testing, in which a specimen is, simplistically, shaken 
in order to impart a high level of energy, can produce a high signal to noise ratio, and can 
induce random excitation, which can average out slight nonlinearities and extract linearized 
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eigenparameter parameters. However, shaker testing is not practically employed in the field 
on relatively large structures, and can be cost prohibitive to conduct. Single impact hammer 
testing addresses the shortfalls of shaker testing, in that it is portable and inexpensive to 
conduct. However, single impact hammer testing falls short where shaker testing is strong, 
in that low energy input of single impact hammer testing produces a low energy input, a low 
signal to noise ratio with no randomization. To address the need for a system which 
combines the advantages of shaker testing and single impact hammer testing, a Random 
Impact Series method for hammer testing is presented which yields a high energy, high 
signal to noise ratio, random system. A novel stochastic model is developed to simulate the 
random impact series produced manually and to generate a random impact series for a 
specially designed random impact device. 

[125] Figure 31 shows a random impact device 1550, according to one embodiment 
of the invention. Random impact device 1550 includes random signal generating unit 1560 
and a random impact actuator 1570 with impact applicator 1580. The impact applicator 

1580 preferably includes a sensor 1581, such as a force transducer, attached at its tip. Sensor 

1581 is preferably configured to send data to the spectrum analyzer in stiffness parameter 
unit 103 in order to obtain mode shape information, as discussed above in connection with 
Fig. 1A. Sensor 1581 can be coupled to the spectrum analyzer in any manner know in the art 
including, but not limited to, wire, optical fiber or even a wireless connection. Random 
impact signal generating unit 1560 generates outputs 1590 and 1591 to random impact 
actuator 1570. It should be appreciated that outputs 1590 and 1591 could be coupled to the 
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random impact actuator 1570 in any manner known in the art including, but not limited, to 
cables, wires, optical fibers, wireless communications, and so forth. 

[126] Output 1590 corresponds to the amplitude \\f h as discussed below. In 
particular, random signal generating unit 1560 outputs a value corresponding to ij/i, as 
discussed below. Output 1591 corresponds to Ti, as discussed below. 

[127] Impact applicator 1580 is shown, for purposes of illustration, as impacting 
structure 1600 in a reciprocating motion. Impact applicator 1580 is shown to have an 
impact path 1610, such that when a structure 1600 lies within the impact region 1610, impact 
applicator 1580 impacts structure 1600 with a force of random amplitude that arrives at a 
random time, as discussed below. The impact applicator 1580 and impact region 1610 as 
shown in Figure 31 is one possible embodiment, and other shapes and impact regions may 
be used while still falling within the scope of the present invention. 

[128] Although random signal generating unit 1560 is shown to output two outputs 
1590 and 1591 that correspond to \\r t and ti below, it should be understood that signal 
generating unit 1590 could output any signal or signals that ultimately results in random 
impact actuator 1570 driving impact applicator 1580 to impact structure 1600 with random 
arrival times x t and random amplitudes y,. 

[129] As discussed above, vibration information is obtained by sensor 110, and is 
sent to stiffness parameter unit 103 via sensor coupler 113. The stiffness parameter unit 103 
sends stiffness parameters to the damage information processor 123. While the random 
impact device can be used for damage detection purposes as discussed above, it can certainly 
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be used just for obtaining the natural frequencies and/or mode shapes of the structure for 
modal testing purposes. 

[130] Experiments conducted on lightning masts by the applicant confirm that 
multiple impact testing performs better than single impact testing when there is wind 
excitation to the masts. Results are shown for a 65 foot tall mast with a 5 foot spike, as 
shown in Figure 32, referred to hereinafter simply as a seventy foot tall mast. The mast has 
two constant cross section schedule 40 pipes of equal length. It can be seen from the results 
shown in Figures 33A - 33B that the multiple impact test has much better coherence, is less 
noisy away from the resonances, and can pick up the modes that were missing in the single 
impact test. The multiple impact test is better at exciting the lower modes that can also be 
excited by the wind, which can be seen by comparing the frequency response functions 
(FRF) between 0 and 30Hz. It can been seen from the FRF for the multiple impact tests 
that there are some modes that are missed at near 4 Hz and 7 Hz with the single impact test, 
and the mode at 10 Hz is much improved compared to the single impact test. 

[131] A stochastic model will now be discussed which describes the random impact 
series Fi - F n modeled as a Poisson process. The Poisson process is one of a general class 
of processes that arise in problems concerning the counting of events in the course of time. 
The force pulses in the series are assumed to have an arbitrary, deterministic shape function 
and random amplitudes and arrival times. The force signal in a finite time interval is shown 
to consist of wide sense stationary and non-stationary parts. The expectations of the average 
power densities associated with the entire force signal and the stationary part of the signal 
are derived and compared, and the power spectral density is related to the average power 
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density and the autocorrelation function associated with the stationary part of the force 
signal. Numerical simulation is conducted to validate analytical predictions, and a 
relationship between the Fourier transform and the discrete Fourier transform used in a 
numerical simulation is produced. Experiments on the four bay space frame, as shown in 
Figure 10, validated the distributions of the random variables and the Poisson process. 

Stochastic Model of a Random Impact Series 

[132] A random impact series is modeled here as a sum of force pulses with the 
same shape and random amplitudes and arrival times: 

*(') = £ V,y (' - h \ ' e (0, oo) (49) 

where t is time, x(t) is the time function of the force signal, r. e (0, t] is the random arrival 
time of the /-th pulse, and yit-Tg) is the deterministic shape function for all the pulses, y/ 1 
is the random variable describing the amplitude of the /-th pulse, and N(i) is the number of 
the pulses that have arrived during the time interval (0,f] and is modeled as a Poisson 
process with stationary increments. All the pulses are assumed to be of width Ar, and 
yit-Ti) satisfies > y(r-r / ) = 0 if t <r i and f>r # .+Ar. 

[133] Since a finite time record is used for the force signal in modal testing, we 
consider only the pulses that arrive during the time interval (0,7] of length T, as shown in 
Figure 34. Note that y(0) and y(Ar) do not have to be equal to zero in this model. 
Because the pulses arriving after time t (t < T) will not affect the force signal at time /, the 
random process N(t) in (49) can be replaced by N(T) . Also, if a pulse arrives at time T, it 
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will vanish at time T + At. The last pulse in Figure 34 arrives at time r N(T) and ends at a 
time between T and T + Ar . To include completely this last possible pulse, we consider the 
time interval t e (0, T + At] . Equation (49) is rewritten as 

N{T) 

*(')= X>,H'- T .)> (50) 

t,. e(0,T] and te(0,T + Ar] 
[134] For the Poisson process N(t) with stationary increments, the probability of the 

, , «-*(*)' . 

event {N(t) = n}, where « is an integer, is P (W) = where * ls the constant 

arrival rate of the pulses. By replacing / with T, the probability of the event {N(T) = n} is 

All the arrival times T i9 where i = 1,2,3,", N(T), are identically distributed, mutually 
independent random variables. Because the arrival rate of the pulses is constant, the arrival 
times r i are uniformly distributed in [0, T] with the probability density function 

T 0< - T< - T (52) 

0, elsewhere 

Similarly, if/ . (i = 1,2,- N(T)) are identically distributed random variables, which are 
mutually independent and independent of the distribution of the arrival times r i . While the 
distribution of ^.(z = l,2,-- 5 N(T)) is not used in the subsequent derivation, it is assumed 
that if/, satisfy the Gaussian distribution with the probability density function 

p (il,) = -—e 2 ° 2 , 0<^<oo (53) 
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where// = E[%] is the mean and a 2 =£[^ 2 ]-^V/l is the variance of y/ r While the 
distribution in (53) is not used in the analytical derivations, it is used in the numerical 
simulation and validated experimentally for a manually applied random impacts on a four 
bay space frame, as shown in Figure 10. 

Force Spectrum and Its Expectation 

[135] The force spectrum is the Fourier transform of the force signal in (50): 

Z M'-w-"* = I « £ y{t-T,y ja, dt (54) 

1=1 J '~1 

where F denotes Fourier transform, *(./'<») is the Fourier transform of x(r), « is the 
angular frequency, and; = ^. Let u=t-r., then t=u+r t and *//=</«. Equation (54) 
becomes 



N(T) 



xuco) = Z v^~ tor/ i ' y W (55) 
i-i 

Since the integral in (55) is a deterministic function, the expectation of the force spectrum is 



E[X(jo>)] = E 



N(T) 

Z w*" 1 

1=1 



(56) 



Using the expression for conditional expectations, | U)] = E(V), where U = N(T) 



N(T) 

and ^ = Z V/ e ~""" r ' ' we have from ( 56 ) 



N(T) 



Z 



= ±P lN) (nJ)E 



N(T) 



-im 



1=1 



£ y(u)e- jm, du 
£ T y(u)e- J °>"du 



(57) 



Since y/ e~ >r ' (1 = 1,2,— ,n) are independent of each other and iff, is independent of e~ jar ' , 
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we have 



± ¥l e~^ 1 = ±E[ ¥i e-^ ] -tE[ Wl ]E[e-»" ] 

_ ,=1 J i-l /=! 

Since ^ (i = 1,2,- ••,«) are identically distributed random variables and so are r f 
(/ = 1,2,- ••,«), we have from (58) 

t^'] = nE[ Wi ]E[e-^] 
Substituting (59) into (57) yields 

n=0 L 

Substituting (51) and (52) into (60) yields 



/j=0 



(AT)"e- XT 



nE[y/ x ] 



j) 



»-i -at 



jco * 

where Taylor expansion of e' 17 has been used. 

Average Power Density of x(t) in [0,T + Ar] and Its Expectation 

[136] The average power density of the force signal in (50) is defined as 
X(jco)X\jco) 



S,{co) = 



T + At 



where X* (jco) is the complex conjugate of X(jco): 



W < r > Mr 
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Substituting (55) and (53) into (62) yields 

*UD mil 



m=I /=1 



T + At 



_ 7 m=l <=1 



(64) 



T + Ar 



(65) 



(66) 



Let k-u-v , then w = k + v and du = dk. The double integral in (64) becomes 

+ ^ r tffvJ[ V+A V(v + ^)^(v> _ya '*<3?A: 
Interchanging the order of integration in (65) yields 

f 1 ^y{u)y{v)e- jm ^dudv= ljkl T k y(v + k)y{ V y j °> k dv 

+ ^ dk^~ k y{v + k)y( V y imk dv 

Let v + k = y , then v = y - k and dy = dv. The first integral on the right-hand side of (66) 
becomes 

J° Ar dk £ r y (v + k)y(vy jak dv = £ r dk ^ z+k y{y)y(j - k)e- Jak dy (67) 
Changing y in (67) back to v and substituting the resulting expression into (66) yields 

f f ' y(u)y(vy M "- v) dudv = £ r f H %(v+ 1 k \)y(y)dve' ja>k dk (68) 

Mr-|A| 

Noting that y(y+ \ k \)y(y)dv is an even function of k, we have 

£ H %(v+ 1 k \)y{v)e j0)k dk = f H %(v+ 1 * |);/(v)cos(<y*)</£ (69) 
Substituting (69) into (68) and substituting the resulting expression into (64) yields 

Si{co)= tt*(& r'^( v+ i*i)^( v )^ c ° sfi> ^) w f zw^-^ (70) 

Since the double integral in (70) is deterministic, we have 

1 /Mr Mr-lAl \ N(T) N(T) 

N(T)N(T) 

Since £ ^>n[;fi>(r, -r m )]} = 0, y/ i (i = 1,2,---,N(T)) ate identically distributed 



m=\ i*=1 



random variables, and so are r. (/ = 1, 2, ■ • • , N(r) ), we have 



N(T)N{T) 
m=1 /'=! 



N(T)N(T) 

£ E^, cos«(r, -r ffl ) 
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= E 



N(T) N(T) N(T) 

£ + £ £ V/ cos °> ( t > <- T n, ) 

1=1 m=l /=! 



(72) 



= E{N(T)tf +[N 2 (T)-N(T)]tf cos^r, -r 2 )} 



Since r, and r 2 are independently, uniformly distributed random variables in [0,7], the 



probability density function of r, - r 2 is 
'T-\t\ 

-T <t <T 



-2 ' 



(73) 



0, 



elsewhere 



Using (51) and (73) in (72) yields 

~N(T)N(T) "I oo ^ 



m=l i=l 



n=0 



+ Z ^} ( r ) ( " 2 - * ) £ [Vi 2 ] £ cos(^r)p n _ r2 (r)dr 



(74) 



Substituting (74) into (71) yields 

^[5 1 (^] = ^[Cr%(v + W)^(v)^vcos(^] 



x\2A>E>M l - C0 ^ T h ATE[tf]} 



(75) 



6; 



Mean and Autocorrelation Functions of x{t) 

[137] The first-order cumulant function of x(t) in (50), *i [*(/)], is equal to its 
mean function, 2s [*(/)]. The second-order cumulant function of x(t), k 2 [x(/ 1 )jc(^ 2 )] , is 



related to its autocorrelation function, x (f 2 )] through 

£ [*h)*fe)] as * 2 [*h)*fe)]+*.[*h)]'r,[xfe)] 



(76) 



where /, and t 2 are any two time instants in [0,T + Ar]. Following the derivations, the first- 
and second-order cumulant functions of x(t) are 
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k x [x(f )] = E[x(t)] = AE[y/ x ] £ y(t-a)da 
K 2 [xMx(t 2 )]nK„(t i ,t 2 ) = ZE[tf)[ y (t l -a)y(t 2 -a)da 

where / € [0,T + Ar] . Let t-a = u in (75), then rfa = -cfa . We have from (77) 



(77) 
(78) 



£[x(0] = A£[^]£ r y(ii)A 



Let 



(79) 
(80) 



(81) 



^(0= [y{u)du 

Noting that y(u) - 0 when u < 0 and « > Ar , we have from (79) 

0<?<Ar 

E[x(t)] = < AE[^]W(At), Ar<t<T 

AE[if/ x ][W{Az)-W(t-T)\ T<t<T + Ar 

where T > At is assumed. Let f,-a = w and t 2 -t l =k in (78), then da = -du and 
/ 2 -or = u + k . We have from (78) 

*« C ,h) = AE [tf ] ||'_ r v( M ) y(u + *)</w (82) 

When \k\>Ar, y(u)y(u + k) = 0 and hence /r„(^ 2 ) = 0 . When 0< A: < Ar , we have from 
(82) for different t, 

AE^y/f^ y(u)y(u + k)du, 0</, <Ar-k 
AE [^]C k y{u)y(u + k)du, Ar-k<t x <T 
AE[tf~\£ T ~ k y(u)y(u + k)du, T <t } <T + Ar-k 

T + Ar-k<t x <T + Ar 
When -Ar < k < 0 , we have from (82) for different t x 

AE[y/f]£ k y(u)y(u + k)du, -k<t } <Ar 
AE[tf~\£ T k y(u)y(u + k)du, Ar<t x <T-k 

AE[tf~\£ T r y(u)y(u + k)du, T-k<t x <T + At 

0, 0<t x <-k 
Let u + k = v in (84), then u = v - k and du = dv. We have from (84) after changing v back to 



(83) 



(84) 
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K xx ('l>0 _ 



ZE[tf]£ +k y(u-k)y(u)du, -k<t x <Ar 



XE [tf ] £ T+k y {u - k) y (u ) du, Ar<t,<T-k 



(85) 

XE\y] ] y(u-k)y(u)du, T-k<t x <T + Ar 

0, 0<t.<-k 
Combining the second equations in (83) and (85), we have for t x and t 2 in [Ar 9 T] 

^ 2 ) = ^E[ W f]^y(u)y(u + \k\)du (86) 
[138] Since by the second equation in Eq. (81), E[x(t)] is a constant for t e [Ar, T], 

and by (86), K„(t x ,t 2 ) is a function of k = t 2 -t x for t x and f 2 in [Ar.7 1 ], x(t) is a wide-sense 

stationary random process in [Ar,T]. Substituting the second equation in (81) and (86) into 

(76) yields the autocorrelation function for t x and t 2 in [Ar,T] 

E[x{ tl )x(h)] = XE[^]^y(u)y(u + \k\)du + X 2 E'[^]W 2 (Ar) (87) ■ 

Fourier Transform of the Mean Function of x(t) and Its Equivalence to E[X(ja>)] 
[139] Applying the Fourier transform to E[x(t)] in Eq. (81) yields 

F{E[x(t)]} = AE[y x ]{ £ r W(t)e- ja "dt + £ W{Ax)e i(i,, dt 

+ £ &T \W(AT)-W(t-T)]e- je "dt) 
The three integrals in (88) are referred to as /, , I 2 , and 7 3 , respectively. Consider first the 
third integral in (88) 

/ 3 = f^ +Ar [^(Ar) - W(t - r)]e- ya "^ (89) 
Let t - T = 6 in (88), then t = T + 0 and dt = d0 . We have from Eq. (89) 

/ 3 = fF(A t)e- j0> e d0 -e imT W{0)e J °" 3 d0 (90) 

Changing 6 in (90) back to / and combining 7, and 7 3 yields 

7, + 7 3 = e- y<Br W(Ar)e- J °"dt + (1 -e"^) W{t)e J °"dt (91) 
Adding 7 2 to (91), simplifying the expression, and substituting it into (88) yields 
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F{E[x(t)]} = XE{ ¥x m-e JoT ) f W{f)e iM dt 

+ e- JaT £ W^r)e- ja, dt + [ t W(*T)e- Jm dt] (92) 

= A^J(l-^)^(A T )[^ + f(^-lK-*] 
[140] We will show that F{£[x(0]} in (92) is equivalent to E[X(jo)] in 

(91). By (80), we have 

0, /e(-oo,0) 

W(t) = - [l{t)dt, /e[0,Ar] (93) 
W{At), fe(Ar,co) 
The integral in (61) can be written as 

y(u)e- jau du = £y(u)du+ ^ y{u){e jou -l)du 

= £ y(u)du - jo) f f v(u) J[ e-^Vvrfw 
Interchanging the order of integration in the double integral in (94) and using (93) yields 
f ' y(u)e im du = £ y(u)du - jco £ dv £ y{u)e j(av du 



(94) 



W(v) . 1 ( 95 > 

JF(Ar) J 

Substituting (95) into (61) yields (91). This shows that the expectation E and the Fourier 
transform F are commutative as both are linear operators. 

r 1 

[141] Equation (92) consists of two parts: the first part, XE[y/^\{\-e Ja )W(At)—, 
is the Fourier transform of the stationary part of E[x(t)] in [At,T], and the second part, 
A£[^,](l-e" >r )^(Ar)J^ r (^^ — l)e~ Ja "dt, is the sum of the Fourier transforms of the 



W{t) 



nonstationary parts of E[x(t)) in [0,Ar] and [T,T + At]. When Ar ->0, since ^7^ _1 is 



finite, — \) e - Je "dt, and consequently the second part of E[X(J<»)], approaches 



W(Ar) 



zero. 



Average Power Density of x(t) in [Ar,T] , Its Expectation, and Power Spectral 
Density 

[142] Since x(t) is stationary in [At J], the average power density of x(t) in 
[At,T] is defined as 
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2V r-Ar 



(96) 



where X s (ja>) = £ x{t)e~ Ja 'dt . Taking the expectation of (96) yields 

Let * = f 2 -f, in (97), then dfc = *ft 2 . Since W 2 )] = ( 97 ) becomes after 

interchanging the order of integration 

4 s ' (•)] - rrs i>' C *- (*)•*** " CL, fi - W( 1 -r^]* M * (98) 

Substituting (87) into (98) yields 



*[*(.)] - «[rf ] CI, P^+Wj-W* 



i- 



:«(*) 


1- 

V 


1*1 


r-Ar 








r-Ar J 



(99) 



Noting that 1 * |) = 0 when | A: |> Ar and %(w + |A:|)>>(w)f/w 

even function of | k \ , we have from (99) 

r / ,n •> i w v 1 -cos <» (r-Ar) 
£[5 2 («)] = 22} E 2 [ ¥ ,]W 2 (Ar)— 



1 — 



1*1 



r-Ar 



is an 



where T > 2 At has been assumed. 



r-Ar 



(100) 



cos cokdk 



[143] The power spectral density of x(t) can be obtained from (100) by increasing T 
to infinity 

S x (a>) = lim S 2 (o>) = 2nX 2 E 2 [y/ x \W 2 {At)S{(o) 



+ AE[tf] £ f _W /( M +|it|)/( M )cos(^)^ 



where we have used 



lim 



l-cos<y(r- At) 



2 sin 



fl>(r-Ar) 



= lim 



llH1 _ U11J — — ^ , (102) 

^ 2 (r-Ar) r->«o nco 2 {T-Ar) 

in which £(□) is the Dirac delta function. The power spectral density can also be obtained 



from (98) by increasing T to infinity 

S x {co) = lim E[S 2 {a>)] = (*)«^<ft =£/U*K^* =^[^ (*)] 



(103) 
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where R^i-k) = R^ik) has been used. Equation (103) is the well-known Wiener-Khintchine 
theorem, which states the power spectral density is the Fourier transform of the 
autocorrelation function. Substituting (87) into (83), and noting that I(u)I(u+ \ k |) = 0 when 
|&|>Ar and the Fourier transform of 1 is 2nS{co) y yields (101). Note that the power 
spectral density is only defined for a wide-sense stationary process with an infinite time 
record. When the mean amplitude of each pulse E[y/ X ] is not equal to zero, there is an 
associated delta function in the power spectral density. 

Comparison of E[S x (a>)\ and E[S 2 (a>)] 

[144] By (100), E[S 2 (co)] consists of two parts: the first part, 

2 2r i 2/ \ l-cos^T 1 — At) 
2X E [y/ x \W (At) — a ) which depends on the arrival rate A, the mean 

amplitude of each pulse E[if/ { ], the total area of the normalized shape function W{At), and 
the time length T-At, describes the average effects of the stationary part of x(t) in [At 9 T] 
and is referred to here as the first-order statistical power density, and the second part, 



* rM*i / / , . ( x |*| " 

v t ~ At ; 



cos akdk , which depends on the mean square 



amplitude of each pulse E[t//f] and the shape function y(Q> in addition to A, T , and At, 
describes the variational effects of the stationary part of x(t) and is referred to here as the 
second-order statistical power density. While E[S x (ca)] in (75) consists of two parts, the first 

Pa "' 2 ^ + ir^ [-Crf r %( V H*IM V )^ VC0S (^*)^*] l ~ co ^ (oT ) f depends also on the 
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shape function y(D) as the shape function is used in calculating the power associated with 
the nonstationary parts of x(t) in [0, At] and [T, 7 + At]. 

[145] When the shape function is a delta function (i.e., Ar -» 0), the nonstationary 
parts of x(t) vanish and £[£,(#)] in (75) can be shown to be equivalent to E[S 2 (o))] in 
(100). By (68) and (69), we have 

JT, f T ^y{ v+ \ k \)y{ v ) dvcosQ)kdk =\^ T y(v)e- Jm dv (104) 

When y(v) = W 0 S(v) , where W 0 is a constant, we have from (104) 



lim ^ y(v + \k\)y(v)dvcosc)kdk = \im^ T W 0 S(v)e- Jav dv 



2 = K (105) 



Since | & |< Ar in (100), we have | k |-> 0 as At -» 0 and hence 

lim 1 1-JAL] = 1 / 106 n 

Substituting (105) and (106) into (100) and noting that W{Ar) = W 0 when y(v) = W 0 S(v), and 
substituting (105) into (75) and noting that 7 + Ar -» T as At -> 0, yields 

£ [5 2 = £ [5, («)] = 2X 2 E 2 ] 1 "°° 2 S 7 f ,r W 2 + A£ [^ 2 ] *F 0 2 (1 07) 

By (104) the power spectral density in (101) can be expressed as 

S x (co) = 2nX 2 E 2 [ ¥x ]W 2 (A t)S{co) + XE [tf ] \Y(jo))\ 2 (1 08) 

where Y(jco) = ^ y(t)e~ Ja "dt . When T -> 00 in (75), we have by using (104) 

S(a>) = lim £[S,(c>)] = 2aA 2 £ 2 [^][y(./fl>)| 2 5(fi>) + A£[> 2 ]|r(»| 2 (109) 
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Equation (109) has a slightly different form from that of (108) because the power associated 
with the nonstationary parts of x(t) is included in (109). When y(t) = W 0 S(t) (108) and (109) 



reduce to 



S x (co) = S(co) = 2ttA 2 E 2 [^ ]W 2 8(co) + AE [y/ 2 ] W 2 



(110) 



Equation (110) can also be obtained from (106) by letting T -> oo and using 

2 sin 2 



2 coT 



]- cos coT ,. 
hm ; = hm 



7 ^ = S(co) 



Examples and Numerical Simulation 



(111) 



[146] When the shape function of the pulses is represented by a half sine wave, i.e., 
^(O = sin^][//(0-i/(/-Ar)] (112) 
where H(0} is the Heaviside function, we obtain by using (60), (64), and (103) 



E[X(jco)] = -j- L 

]CO\n -co At J 



(113) 
(114) 



(115) 



£[^,H] = 2 f A ^ 1+ 2 C ^ A ^ ATE^)- 2XE ^ C 2 OSG)T -^ 
1 IV n (co 2 At 2 -7t 2 ) 2 (T + At)[ ri ' co 2 

r , xn , f"-27 , ^- 4 -2rcos((yAr)^ 4 +(4cos(6jAr)^ 4 +^ 4 )Arl 

[+(8sin(<aAr)a»r 2 + 2rcos(<yAr)<yV + 27VV)Ar 2 ] 
(<y 2 Ar 2 -;r 2 ) 3 (r-Ar) 
| [(ft; 4 Ar 5 -4cos( 6 ;Ar)ft;V-2a;V)Ar 3 ] ^ 8^ 2 (^)Ar 2 q-cosQ>(r-Ar)) 
(<y 2 Ar 2 -;r 2 ) 3 (r-Ar) l+ ^V(r-Ar) 



Consider next the normalized shape function y(t) shown in Figure 35 with unit maximum 
amplitude. It is obtained by averaging a series of normalized force pulses from impact tests 
on the four-bay space frame as shown in Figure 10. There are 21 sample points in the shape 
function, which are connected, as shown in Figure 34. Other parameters used are T = 8 s , 
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Ar = 20x7/1024 = 0.15625 s where h = 771024 = 0.0078125/5 is the sampling interval, 
/l = 4.14/s, E[^ ] = 0.8239 N, and E[y/ 2 ] = 0.7163 N 2 . The curve for E[x(t)] in the time 
interval from 0 to 8.15625 s, shown as a solid line in Figure 36, is calculated using (81). It is 
seen that E[x(t)] increases from 0 to 0.0352 Nin the first 0.15625 s, remains at 0.0352 N 
from 0.15625 s to 8 j-, and decreases from 0.0352 N to 0 in the last 0.15625 s. The curve for 
20 log | E[X(jo))] | in the frequency range from 0 to 50 H%, shown as a dotted line in Figure 
37, is calculated using (61). The curve for 101og{£[S, (y'<y)]} in the same frequency range, 
shown as a solid line in Figure 38, is calculated using (75) . The curve for 10 log {£[.!>, (y'<y)]} 
with A = \/s and the other parameters unchanged is shown as a dashed line in Figure 38. It 
is seen that £[5, (y'ty)] increases by 4.14 to 15.6247 times in the frequency range shown 
when X is increased from ^ =1/5 to \ =4.14/5. This result can be shown by using (75) 

W] U _ tfE 2 [^]T 2 + A,TE[tf] ^E 2 ty x ]T + \E[tf] _^ ^ 
\\mE{S,(jo>)]\ x=h X^E'W^+X.TEW 2 } V*V,F + <W, 2 ] 

iszzz ^L = A= 4 . 14 (117) 

lim£[S, (;«»)] 4 

This shows that a larger arrival rate X would increase the energy input to the structure over 
the entire frequency domain. 

[147] Numerical simulation is undertaken next to validate the analytical predictions. 
The random number N(T) satisfying the Poisson distribution in (51) with A = 4A4/s and 
T = 8 s is generated using MATLAB . Similarly, the random numbers corresponding to the 
random variables r i (i = l, 2, N{T)) y satisfying the uniform distribution in (52), and the 
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random numbers corresponding to the random variables ^ (i = l, 2, — , N(T)), satisfying 
the Gaussian distribution in (53) with // = 0.8239 N and 
a 1 = 0.7163 -0.8239 2 N 2 = 0.0375 N 2 , are generated. Using the shape function constructed 
earlier, a sample function of x(t) in (50) at time t = rh, where 

r = 0, 1, • • , R - 1 ( R = r + AT = 1 044) , denoted by x r , can be obtained. The discrete Fourier 

h 

transform (DFT) of the time series {x r } is calculated using MATLAB. The DFT of the 
series {x r } is defined by 

i r-i ,i*v 

*,45> (118) 

where q = 0,1, • Equation (118) is an approximate formula for calculating the 

coefficients of the Fourier series of a periodic function whose values in the period 
[0, T + Ar] are given by those of x{t) : 

1 rr+Ar -7 



X =-J—r T x(t)i J ^dt (119) 
9 J + Ar* 



The Fourier components AT ? correspond to harmonics of frequency <o q = • Recall 

that the Fourier transform of x(t) in (50) is given by 

X{ja>)=l^x{t)e j °>dt (120) 

Let co = <y = g in (120) and compare the resulting expression with (119), we find that 
' T + Ar 

* in (118) multiplied by T + Ar provides an approximate value of X{j<o) at frequency co q . 
Similarly, X q X\ multiplied by T + Ar provides an approximate value of 

s ( W \ = X U°>) X 0'^) at frequency <y fl . By averaging 5000 sample series of {x r } , we obtain 
1 T + At 
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the curve for E[x(t)], shown as a dotted line in Figure 36. By averaging 1000 sample series 
of {201og[(r + Ar) | X q |]} , we obtain the curve for 20 log | E[X(jco)) \ , shown as a solid line 
in Figure 37. By averaging 100 sample series of {101og[(r + Ar)|* 9 ** |]}, we obtain the 
curve for 10 log {£[£,(>>)]}, shown as a dotted line in Figure 38. The numerical results are 

in good agreement with the analytical ones. 

The stochastic model was experimentally validated for an experimenter conducting 
manually a random series of impacts on the four bay space frame as shown in Figure 10. 
One hundred ensemble averages were used. The experimental probability density functions 
of the arrival time, the number of arrived pulses, and the pulse amplitudes are in good 
agreement with the analytical values, as shown in Figures 39-41. 

[148] Thus, the system and method for detecting structural damage and the random 
impact series method as embodied and broadly described herein can be applied to an 
unlimited number and type of structures to provide automated, reliable damage detection 
and assessment and to conduct modal testing. This system could be further automated to 
conduct periodic tests and provide results to a centralized monitoring section. Regular 
health monitoring of these types of structures could provide additional protection against 
potential failure, as well as a characterization of usage and wear over time in particular 
environmental conditions for predicting useful service life. 

[149] The foregoing embodiments and advantages are merely exemplary and are not 
to be construed as limiting the present invention. The present teaching can be readily 
applied to other types of systems. The description of the present invention is intended to be 
illustrative, and not to limit the scope of the claims. Many alternatives, modifications, and 
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variations will be apparent to those skilled in the art. In the claims, means-plus-function 
clauses are intended to cover the structures described herein as performing the recited 
function and not only structural equivalents but also equivalent structures. 

i 
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